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This report describes < computer progrem developed at the Douglas Air- 
erait Division of the Mc Donnell Douglas Corporation, Long Beach, Cali- 
fornia, The development of the Douglas Hypersonic Arbitrary~Body Aero- 
dynamic Computer Program was started in 1964 and greatly expanded in 
aubsequent years under sponsorship of the Douglas Independent Research 
and Development Program {IRAD}. From August 1966 to May 1967 the 
program development was continued under Air Force Contract No. F 3361567 
1603. This contract was administered under the directien of the Aerc- 
nautical Systems Division, Directorate of Analysis, Wright-Patterson Air 
force Base, Ohio by Mr. R. K. Milis, Project Engineer (ASBED-30}. 

The product of this work was the Mark Il version of the prédgram as released 
for use by government agencies in May 1967, The latest version of the 
program as presented in this report (the Mark Ui version} is an extensively 
revised version of the earlier Mark If program, This version has been 
prepared az a result of both 1967-68 Douglas IRAD work and another Air 
Force contract {F336i5 67 C 1602). This contract was administered by the 
Air Force Flight Dynamics Laboratory, Flight Mechanics Division, Gas 
Dynamics Branch, Mr. Valentine Dahlem, Praject Engineer (FDMG). 


At the Deuglas Aircraft Division this work was conducted under the 
direction of Mr. A. E. Gentry as Principal Investigator. A number of 
peopie contributed to the various phases of this work for which the author 
ig grateful. Mr. BD. N. Smyth provided valuable consulting services in 
many phases of this work and prepared the new skin friction techniques 
incorporated in the Mark Ii version. Mr. W. R. Oliver's work in applying 
the various versions of this program to practical design problems con- 
tributed both in program design and in program validation. Others partici- 
pating in this work include Messrs. G. D. Buell, J. L. pee N. F, 
Wasgon, and B. G, Wilson. 


Special appreciation is extended to the various users of the earlier versions 
of this program for their valuable suggestions in a number of areas and for 
their efforts in adapting and running earlier versions of the trogram on the 
different types of computers. These include Messrs. Fred White, Jr. (Afr 
Force ASBED-30), Don Shereda (Air Force FDMG}, Ralph Carrnichael and 
Charies Castolano (NASA Ames), C. L. W. Edwards (NASA Langley), 
Ralph Grahm (NASA Houston), Ray E. Aley (Lockheed Electronics Co., 
Houston), and R. E. Finch, A. W. Marziane, and J. H. Kainer (Aerospace 
Corp. ). 


This computer program and documentation report were released for generai 
use by the author and by the Guidance and Control Section, ASBED-30, 
Wright-Patterson Air Force Base, in April 1968, This program and report 
are provided in the interest of information exchange. Responsibility for 

the contents rests with the author or organization that prepared it. 


The dis {bution of computer program decks for the Mark III version is 
handied “he author. 
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ABSTRACT 


This report describes a digital eomures program system that is cagabie 
of calculating the hypergonic aerodynamic characteristics of complex 
three-dimensional shapes, The outstanding features of this program are 
its flexibility in covering 4 very wide variety of prablems and the multitude 
of program options available. The prugram ie a combination of techniques 
and capabilities necessary in performing a complete aerodynamic analysis 
of hypersonic shapes, These include vehicle geometry generation and 
description, visual graphics necessary in handling geometry data and in 
preparing plots of the final aercdynamic data, aerodynamic calculations 

of surface pressures and skin friction forces, and the integration of these 
forces to give ali aerodynamic coefficients and stability derivatives. , 


The geometric description techniques in this program provide the capability 
of handling completely arbitrary three-dimensional shapes. The procedure 
developed to check the accuracy of the geometric data uses @ computer and 
automatic recorder to draw pictures of the vehicle viewed from any angie, 


The pressure calculation methods provided within the program inciude 
modified Newtonian, blunt-body Newtonian-Prandti- -Meyer, tangent-wedge, 
tangent-cone, shock-expansion, Prandtl-Meyer expansion, blast wave, 
modified tangent-cone, boundary-layer induced pressures, free-molecular 
flow, and a number of empirical relationships. The pressure calculation 
method most suitable for each component of the vehicle is specified by the 
aérodynamiicist. Viscous ferces are also calculated and include viscous-~- 
inviscid interaction effects. Skin friction options include the Reference 
Temperature and the Reference Enthalpy methods (for both laminar and 
turbulent flow), the Spalding-Chi method (turbulent), and a special blunt 
body skin friction method. Control surface deflection pressures, including 
separation effects that may be caused by thm deflected surface, are also 
calculated. 


The program has been used to study a wide variety of hypersonic vehicle 
shapes including hypersonic cruise aircraft. air-breathing booster air- 
craft, blunt lifting reentry bodies, high L/L) reentry vehicles, blunt 
reentry capsules, rocket boosters, reentry warheads, and satellite shapes. 


The program is documented in two volumes. Volume I is primarily a 


User's Manual, and Volume II contains the Program Formulation and 
Listings. 
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SECTION I 


INTRODUCTION 


The objectives of the research work that led to this program were to (1) de- 
velop methods for determining the aerodynamic force characteristics of 
hypersonic vehicles regardless of the vehicle shape or flignt condition, 

(2) program those techniques that required digital computer capability for 
practical and efficient application, and (3) verify these techniques by com- 
paring the analytical results with test data. 


At the start of this research project a list of guiding objectives was estab- 
lished to insure successful completion of the work. Major features desired 
in the final analysis system would: 


ie Provide the ability to analyze completely arbitrary 
three-dimensional shapes. 


2. Provide a component build-up capability where each 
vehicle component may be of arbitrary shape. 


3s Include a number of force analysis methods so that 
the system would have the widest possible applica- 
tion to various vehicle shapes and flight conditions. 


4, Provide the capability to use the best force calculation 
method for each vehicle component. 


5. Provide methods for analyzing simple shapes within a 
minimum time period. 


6. Develop a total analysis system framework that is 
adaptable to continued improvement and expansion. 


The initial phase of this work was started in late 1964 and continued in 1965 
as part of a Douglas Independent Research and Development Study. During 
that time a general arbitrary body force analysis approach was derived for 
hypersonic vehicles, the important basic components of the computer sys- 
tem were written and checked out, and the system demonstrated by applica- 
tion to several vehicles of completely arbitrary shape. Alli of this was ac- 
complished under a very modest work effort. 


During 1966 this work effort was expanded slightly and new capability added 
to the program system. This included the incorporation of several] new 
force calculation methods, and the expansion of the force program to cal- 
culate vehicle static and dynamic stability derivatives. 


In August of 1966, the program development was continued under Air Force 
Contract F33615 67 C 1008. This work included further expansion ot the 
force calculation methods, the addition of new geometry description fea- 
tures, the incorporation of control surface derivative calculations, the 
consolidation of all the system components to form one large program, and 
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the preparation of complete program documentation information, The final : 
program resulting from this work was identified as the Mark I version of ar 
the Hypersonic Arbitrary- Body Aerodynamic Computer Program. 


During 1967 -1968 the program development was continued under Douglas 
IRAD and another Air Force Contract (F3361§ 67 © 1602), During this 
period a number of major program modifications were accomplished in- 
cluding the addition of a number of new pressure calculation options, ex- 
tensive revision to the skin friction parts of the program, and conversion 
of the program to operate on different computers. The program resulting 
from this work is identified as the Mark liI version. 


Throughout this report it will be assumed that the reader is familiar with 
the contents of Volume I, the User's Manuai. Discussions of earlier 
versions of this program are given in References 1 and 2. 


Both Volu:nes land II of this report are essentially revised editions of 
the earlier Mark Il program report (Reference 2). The differences be- 
tween the Mark III and Mark Ul reports reflect the modifications and new 
capabilities provided by the latest version. 


This report contains descriptions of the analysis techniques used within 
the program. Throughout these discussions an attempt has been made to 
maintain mathematical notations cousistent with the appropriate refer- 
ences involved, This will assist the reader in comparing the approaches 
with the original reference material at some slight loss in continuity 
within the present report, This policy has also been used in the selection 
of many of the program variable names. 


The program source language itself has been written in a manner so 

that the general flow of the program logic is easy to follow, even by the 
reader unfamiliar with FORTRAN. This has been accomplished by a 

very liberal use cf comment statements throughout the program. This 
approach, together with the availability of machine-produced program 
flow charts, and the general widespread knowledge of the hasic FORTRAN 
language, makes it unnecessary to include detailed equation-by-equation 
descriptions of the program content. Instead, it is only necessary to give 
a general mathematical description of the approach used. 


Appendix A to this report contains the source language listings of the 
program, machine-produced flow charts,-and the definitions of the pro- 
gram symbols. Appendix B contains one complete alphabetical listing 
of all the prograra variables and their definitions. Appendix ~ contains 
a description of the program subscripted variable arrays. 
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SECTION II 
PROGRAM FORMULATION 


PROBLEM FORMULATION 


The problem of estimating the aerodynamic characteristics of 


.arbitrary three-dimensional shapes at hypersonic speeds has sev- 


eral salient features. The first problem is to construct an accurate 
description of the vehicle geometry, as arbitrary and cornplex as it 
may be, This difficult geometry problem, together with the extremely 
wide range of flight conditions, dictates that many different force pre- 
diction techniques be available for use. The various approaches used 
in calculating the aerodynamic forces on three-dimensional shapes 
differ in that different methods are used to attack these two basic 
problems — the geometry representation problem, and the force cal- 
culation technique problem. 


The most ambitious approach presently being attempted uses the 
method of characteristics. Some of this work is discussed in Ref- 
erences 3 and 4. Undoubtedly, this is the eventual ideal approach 
to the calculation of forces on three-dimensional shapes at high 
speeds, However, present mathematical techniques and digital 
computer size and speed capabilities prevent application to typical 
preliminary design problems. Present applications of the method 
of characteristics must be reserved for simple bodies of revolution 
at zero angle of attack or important detail design applications where 
large computer times of several hours might be acceptable. 


Other detailed gas dynamics approaches have been used ior very 
simple shapes such as blunt-nosed bodies of revolution (see Ref- 
erences 5 and 6). While this work has considerable value as an aid 
to understanding the chemical and gas dynamics problems associated 
with simple blunt shapes, its application to the general three- 
dimension .«l shape is restrictive. 


Other special techniques, some theoretical and some empirical in 
nature, have been worked out for simple shapes. Notable examples 
are the work for hypersonic flow conditions on delta wings by 
Creager (Reference 7), and McLaughlin (Reference 8}. 


Another interesting approach to understanding the nature of high- 
speed flow is to select special shapes that are arnenable to exact 
theory. In Reference 9 a class of delta wings is considered that 
permits solution using the exact shock wave theory. This work is 
expanded further in Reference !0where upper surfaces are it.cluded 
and complete configurations suggested that are derived from simple 
shock waves and expansions. 


Reference 1] contains an approach that uses parts of conical shock 
waves, This work is also discussed in Reference 12. A typical 
general approach to these and subsequent similar studies may be 
described as follows. Several exact methods have been derived for 


the solution of Now fields about sight circeiur cones at zero angle of 
attack. The shock wave system forrned by this cone is also of a 
conical shape and the streamlines and flow properties of the air 
pussing through this sheck structurs may he readily calculated (Ref- 
erence 12}, The fixst step in deriving x configuration is to define a 
wing leading edgs line that Hes on the surface of the shock wave cone. 
if we then follow the streamlines formed by the original cone frorn 
this leading edge line downstream, we find that we have described a 
surface along which svluficas are already available from the exact 
eone flow calculations, Such an approach has been aged in designing 
and evaluating same high speed airbreathiny engine inlet designs. 
Although this approach again adds to the "storehouse" of knowledge 
on Ingh speed aerodynamics, it certainly doea not give much flexibil- 
ity sither in the design or the analysis pzecess. 


Amore practical avpreach te the vehicle analysis problem was pre-~ 
sented by Hankey in Reference 14. Im this vork it was recognized 
that each general type of vshicle shage {i.e., leading edge, flat 
suxface, cone} required a difisrent methed for accurate force pre- 
Sctien. In this approach, seéiutions were obtAdined by restricting 
the vehicie to a few set combinations of simple geometrical shapes 
and by deriving theerstical and empirical force calculation methods 
for each shape component. The approach gave good results for the 
ciags of vehicles considered {the Dyna-Soar vehecile shape) but was 
limited in its application because of the restricted set of shapes 
availeble. The major contribution, however, was to illustrate that 
the method of force calculation must be tailored to the vehicle 
component shape and to the specific Hight condition involved. 


A similar approach has been in widespread use in both government 
and industry areas. The general procedure is to break the vehicle 
into a number ef small simple shapes and to estimate the aerodynamic 
{crces on each component by available methode. The widest and most 
frequently used force methed is Newtcnian theory. Several reports 
have been published that give results of Newtonian calculations for 
simple geometric shapss. References 15 and 16 contain design 
charts that are useful in evaluating certain geometric components 
(hemispheres, cylindrical leading edges, cones, cone frustums and 
flat plates) by Newtonian theery. The difficulty of this approach is 
that only simple shapes may be analyzed and only elementary force 
analysis method results such as given by Newtonian theory are avail- 
able. For complex shapes considerable time and effort is required 

to define the vehicle in terms of required simple component shapes. 


Reference 17 uses a semigraphical method for solving the arbitrary 
body geometry problem and uses Newtonian theory for the force cal- 
culations. The work by Dahlem jcef ASD) modifies this procedure to 
allow computation on a digital computer (Reference if}. The geomet- 
ric input data consist of coordinates of points on the surface, two 
angles indicating the surface siope and an angle indicating the ovients~ 
tion of the surface point. The inherent difficuity of preparing these 
input data for complex configurations from conventional engineering 
drawings is obvious. 


.at 


A different approach to the geometry problem iz prevented by 

Van Tassali {ef AVCC} in Reference 19, In this approach the Yady 
is approximated by a series of trianglee on the body surface ag - 
eetermined from input coordinate points. In thig case ths force 
calculation methods include Newtonian theary and free-meolecular © 

Siew (completely diffuse reflection), The use of surface coordinate 
points to furm triangular elements still presente 4 difficult probiem 

in, cbtaining accurate and cotrect input date. Li “- 

A more elaborate technique cf surface description has been worked cut” 
at MIF by Coons (Reference 20). In this method the surface is divided 
into 2 number of lasge patches and then each pateh is mapped into a two~, 
dimensional surface and described by a mathematical surfsce fit ; 
technique, The boundary curves of a patch are described by third order 
polynormials requiring only the coordinates and aurface slopee at the 

‘sur corner points in the transferred coordinate system. In a practical 
application the corner poini surface slopes must be calculated by the 
computer since the preblem of measuring them directly from engineering 
drawings would restrict tha usefulnese of this methad. This may be 
accomplished by using 4 nurmber of surface points along the boundaties 
to more commiletely describe the shape. After the boundary curves are 
defined in the transiormed systern the covrdinates of a peint on the 
surface may be calculated by using special biending relationships that 
serve as weighting functions to properiy account for the influence of each 
boundary. Reference Zl contains a brief desription of a similar 
technique used ky the Boeing Company. 


In most of the referencee cited above simple Newtcnian theory was 
aged te calculate the serndynamic forces, While this method may 

give acceptable resuits for some shapes at very high Mach numbers, 
other matheds must be available before a truly arbitrary body analysis 
system is available. An equally importart problem is the transcrip- 
tion of the geometric description of the vehicle from engineering draw~ 
ings and the preparation of this information in a form acceptabie tc the 
digital computer. Regardiess of what surface description method is 
used, checking of the volurninous and often complex input data poses 

a difficult problem. Since a siugle vehicle may he cornposed of 
several different clagses cf general shapes (such as a blunt nose 

body, together with a thin wing) the arbitrary body anaiysis system 
must have the capability of using different force analysis techniques 
and irequently a different grid work of input surface data on each 
component of the vehicle. 


For hypersonic speeds the most frequently used force analysis method 
is Newtonan theory. In this theory the pressure on each surface point 
of the vehicle is only a function of the angle that the surface makes with 
the free stream flow. Several other methods for calculating inviscid 
pressure also depend only on the local inclination of the surface. These 
include tanygent-wedge and tangent-cone, and several other empirical 
methods worked out for special uses such as Newtonian plus Prandtl- 
Meyer expansion. Another important m2thod for use at high speeds is 
hypersonic shock-expansion. 
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With the foregoing discussion as background informaticn we will now 

_ proceed with the approach used to solve the general probiem that is 
the subject of this research project — the development of methods for 
calcuisting the hypersanic, aerodynamic characteristics of completely 
arbitrary three-dimensional shapes. 


In the following pages we will first present a brief introduction vo the 
yrogram organization. This will be followed by a more detailed deriva- 
tion of the analygia methods used in each component of the system. The 
computer program system described in the following pages is capable 
of calculating the aerodyaamic forces on completely arbitrary three- 
dirnensional shapes at high supersonic and hypersonic speeds. 


All pragram components ef this system are written in FORTRAN. 


PROGRAM COMPONENTS 


Thies program is written on a completely modular basis to facilitate 


checkent and modification activities. The Mark Ill pregram contains 
five major program components: the Aerodynamic Program, the 
Picture Drawing Program, the Output Data Plotter Program, the 


Auxiliary Geometry Generation Program (Slab Delta Program), and 
the card punch routine. . 


Jn the early develcprnent phase each of these components was actually 
a completely separate computer program. The programs could, how- 
ever, be run back-to-back with the output from one program saved on 
tape for use by the next. In the Mark lilprograrn all of these compo- 
nents are combined under the control of a small executive main pro- 
gram. However, to simplify the discussions in this report, each of 
these components is still referred to as a program even though it is 
really just a subroutine in the overall program. 


The general activities carried out by the Arbitrary-Body Program 
System are involved with one of three basic tasks: (1) the preparation 
of geometry data, (2) the calculation of aerodynamic characteristics, 
and (3) the preparation of graphic output data. The computations per- 
formed in support of each of these system tasks are discussed in the 
following sections of this report, 


GECGMETRY 


This program contains several different methods for describing three- 
dimensional geometric shapes. These methods provide the {flexibility 
required to analyze a variety of shapes ranging from very simple sur- 


faces to the most complex forms. The program geometry options provided 


are (1) the surface element method, (2) the elliptical surface generation 
method, (3) the parametric cubic method, and (4) the slab delta geometry 
generation method. If desired, all of these methods could be used in 
describing a single «.ehicle shape. 


The principles involved in the application of each of these geometry 


methods are discussed in detail in the User's Manual and need not be 
covered here. The principal mathematical techniques, however, are 
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important from the programming standpoint and will be discussed on 
the following pages. 


The Surface Element Geometry Method 


The basic geometry method used by this program is the surface 
element or quadrilateral method. This method was developed by 

J, L. Hess and A. M. O. Smith for the Douglas Three-Dimensional 
Potential Flow Program (Reference 22). For completeness, certain 
parts of this report will be included in the following discussions, 


The coordinate system used for this analysis is a right-handed 
Cartesian system a3 shown in the figure below. 


SIAGONAL VECTORS T, and 7, 
Tip = Xs-Xr Tiys¥s- i) Tare 2Zs- 2) 
Tar = Xa~-R2 Try = Va-¥2  Ta25 Za~ 22 
UAIT RORMAL heT,xT, 
Ry = TayTiz - TryTaz fh, 2 N/M 
Ny = TisTaz - TaxTte ty = Ry/N 
Nz = TagTry ~ TrxTay fz = HK 


Ne vige hy + Ne 
AVERAGE POINT 
X = (Xy + Xz + Xs + Xa) 
y = 4(Y, + Yoo Y3 4 Ya) 
22 4(Z; + Z2 + Z3 + Z4) 
CORNER POINT PROJECTION DISTANCE 
dy = nfK - Xy) + OV ~ Va) + 2 - 25) 
k=1,2,3,4 
CORNER POINT COORDINATES 
X = X, + nidy ‘ 
Yi = Yu + Myth 
Zi, = Zy + 120_ 





‘ 


¢ 
In the conventional use of this program the vehicle is usually positioned 
with its nose at the coordinate system origin and with the length of the 
body stretching in the negative X direction. The slight inconvenience of 
this negative sign on the body stations has been accepted so that the 
geometric data will be compatible with the Douglas potential flow pro- 
gram {Neumann Program). 


The body surface is represented by a set of points in space. These 
points are selected on the body surface and are usec by the method to 
obtain an approximation to this surface that is used in subsequent cal- 
culations. If the four related points of each set are connected by 
straight lines we may obtain a picture of how the input surface points 
are organized to descrihe a given shape. This has been done in Fig~ 
are 1, The input scheme has been designed so that each point need 
only be input once even thougn it may be a member cf as many as 
four adjacent sets of points. This is accomplished by the use of an 
additicnal parameter for each point besides the X, Y, and Z values. 
This parameter (known as the status flag) indicates whether a point is 
a continuation of a column of points (STATUS = 0), the beginning of a 
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Figure 1, Output from Perspective Drawing Program 
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new column of points (=1), the first point of a new section of elements 
(=2), or the last point input for the shape {=3}. 


Ae may be seen from the drawings made by the Picture Drawing Pro- 
gram, the different areas of a vehicle may require a different organiza- 
tion and spacing of surface points for accurate representation, Each 
such area or organization of elements is called a section and each sec-~ 
tion is independent of all other sections. The division of a vehicle into 
a given set of sections may also be influenced by another consideration 
since the force calculation program may be made to calculate the force 
contributions of each section separately, using different calculation 
methods. 


The input surface points are not sufficient in themselves for the force 
calculations. Each set of four related points which form an individual 
element must be converted into quantities useful to the program, This 

is accomplished by approximating each element area of the vehicle by 

a plane quadrilateral surface. Since we are using four surface points 

to form an element, no single surface will contain the points themselves. 
Also, adjacent plane quadrilateral edges will not necessarily se co- 
incident. With a sufficiently small size of the surface elements this will 
be of no consequence in the end results. 


The mathematical technique used in converting an input set of four 
points into a plane quadrilateral element is described below. The 
figure below gives a representation of the input element points with 
each point identified consecutively around the element by the sub- 
scripts |, 2, 3, and 4, respectively. 
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We may now obtain a new vector N (and its components) by taking 
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identifies the coordinates as input coordinates. 
next form the two diagonal vectors Tj and Tz. The components of 
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The unit normal vector, A, to the plane of the element is taken as WN 
divided by its own length N {direction cosines of outward unit normal). 


3 
. N 
Zh@ 2h zh 


where _ 
N= nN? +N %4+N2 
; x ¥ z 


The plane of the elernent is now completely determined if a point in 


a 


e this plane is specified. This point is taken as the point whose co- 
& ordinates, X, y, Z are the averages of the coordinates of the four 
; input points. 

: wie,” = 1 | i i i i 

b x = Z xy + X5 + 3 + x4 

i es L (i i 

Y, = tivity t yg t v4 

£ = m i i i i i 

: Z = a [tt ht el 


Now the input points will be projected into the plane of the element 
along the normal vector. The resulting points are the corner points 
of the quadrilateral element. The signed distance of the k-th input 
points (k= 1, 2, 3, 4) from the planeis 


a 


i ares | = i a 
k n(% ~ x) + nF - y,) + ni (2 ~ Zz) k=l, 2, 3, 4 


It turns out that, due to the way in which the plane was generated from 
the input points, all the q,,’s have the same magnitude, those for points 


| and 3 having one sign and those for points 2 and 4 having the opposite 
sign. Symbolically, 


qd). = (- 1) a, k=], 2, 3, 4 
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The magnitude of the common projection distance is called d, i.e., 
a = [a] 


The coordinates of the corner points in the reference coordinate system 
are given by : 


i 
x, = x, + a 
t 2 
Ye = Ve F BLA ks1, 2, 3, 4 
1 
Zr = zy + n, dy. 


Now the element coordinate system must be constructed, This 
requires the components of three mutually perpendicular unit 
vectors, one of which points along each of the coordinate axes of 

the system, and also the coordinates of the ovigin of the coordinate 
system. All these quantities must be given in terms of the reference 
coordinate system. The unit normal vector is taken as one of the 
unit vectors, so two perpendicular unit vectors in the plane of the 
element are needed. Denote these unit vectors ty and t. The 


vector t, is taken as T, divided by its own length Ty: i.@., 


B 


= Ix 

hx i T; 
T 
t = _ly 
ly T, 
t = The 
lz T) 
where 

= 2 - 2 2 

T, = Thx + Ty + Tie 
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The vector tf, is defined by S a nx te so that ite cormponents are 


fax * as - ahy 
toy “ Meee 7 ie 


The vector t, is the unit vector parallel to the x or € axis of the 


element coordinate system, while t, is parallel to the y or 4 axis, 


and nis parallel to the z oz $ axis of this coordinate system. 


To transform the coordinates of points and the components of vectcrs 
between the reference coordinate system and the element coordinate 
system, the transformation matri: is required. The elements of this 
matrix are the components of the three basic unit vectors, ty to» and 
i. Teo raake the notation uniform define 


By Ee Bg Sy Sige 
ny mR YO a 
| A a Yc 233.7 Rg 


The transformation 1». ‘x is thus the array 


ayy #12 213 
43) 392 433 
435 432 433 


To transfer the coordinates of points from one system to the other, 
the coordinates of the origin of the element coordinate system in the 
reference coordinate system are required. Tet these be denoted x , 


Yor 25: Then if a point has coordinates x', y', z' in the reference 


coordinate system and coordinates x, y, z in the elernent coordinate 
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system, the transformation from the reference to the element system 


is 
= ' + ta + = 
x By 1 (x x_) a,oly y_) a. 3(z! z} 


= to = = 
y = ao 1 (x x,) + azoly' pe) + 53 (z! 2.) 


= 23 oe to 
z a3, (x x) + azoly Yo) + agalz z,) 
while the transformation from the element to the reference system is 


' a ; * ee 
x x, + ay 4% t Zo1¥ + ay )% 


y = v + 842% + Bo0V + Big a% 


! ~ 
Zz = z+ ay 3% + ao3Y + 333% 


The corner points are now transformed into the element coordinate 
system based on the average point as origin. These points have co- 
ordinates xy he? Za in the reference coordinate system. Their co- 


ordinates in the elernent coordinate system with this origin are de- 
noted by ou Te 0. Because they lie in the plane of the element, 


they have a zero z or € coordinate in the element coordinate system, 
Also, because the vector th, which defines the x or € axis of the 
element coordinat2 system, is a multiple of the "diagonal" vector 


a 
from point 1 to point 3, the coordinate un and the coordinate N3 


are equal. This is illustrated in the figure below. Using the 
above transformation these coordinates are explicitly 
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: aegorentroid., Se 


(E35, 3) 


¢* = a,,(xt - X)+a,,(y! - J) + a,.(z) - 2) 
k 11k TAO 7 Y 13'% 
K=1, 2, 3, 4 
n* = a,,(x! - x) + a,.,(x! - y) + a,.(z' - 2} 
k 21% 22% 7 ¥ 23'*k 
These 


corner points are taken as the corners of a plane quadrilateral. 


The origin of the element coordinate system is now transferred to the 
centroid of the area of the quadrilateral. With the average point as 


origin the coordinates of the centroid in the element coordinate sys- 
tem are: 


_ J = j * % % % can * 
bo = 3 n* = y# le (1, = nN) + 25 (14 *) 
2 4 
: 1. 
No ia ~ 3° 
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These are subtracted frorn the coordinates of the corner points in the 
element coordinate system based on the average point ag origin to 
obtain the coordinates of the corner points in the element coordinate 
system based on the centroid as origin. Accordingly, these latter 
coordinates are 


b. ‘i B78, 


a * 
Hee FC Mg ale 


Since the centroid is to be used as the origin of the element coordinate 
system, its coordinates in the reference coordinate systen: are required 
for use with the transformation matrix. These coordinates are 


x, = Baz, 65442) 1, 
Yo e Hs 242 bo f ary ‘lo 
Zz, = tay, & ta.3% 


Since in all subsequent transformations between the reference coordinate 
system and the element coordinate system the centroid is used as origin 
of the latter, its coordinates are denoted Xo? Vor %: The coordinates 


of the average point are no longer needed. The change in origin of the 
element coordinate sysiem, of course, has no effect on the coordinates 
of the corner points in the reference coordinate systern. 


The lengths of the two diagonals of the quadrilateral, t] and t>, are com- 
puted from 


t= (bs - &)? 
to” = (4 - &o)* + (ng - 15)2 


The larger of these is selected and designated the maximum diagona! t. 


The body surface area and enclosed volume are determined by summing 
up the contributions of each element. In terms of the coordinates of the 
corner points, the area of the quadrilateral is 


1 : 
A= 5 (€3- &)) (12-74) 


The incremental volume is given by the volume of the parallepepiped 
formed by the element and its projection onto the x-z plane (the x-y 
or y-2z planes would have served equally well). 
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The foragcing procedure may be briefly summarized as follows: 


Each set of four puints is converted into a plane-quadrilateral element 
by the procedure shown in the sketch on page 7. The sormai to the 
quadrilateral is taken as the crosz2 product cf two diagonal vectors 
farmed between opposite element points. The urder of the input points 
ang the manner of 6efining the diagonal vectors is used to ensure that 
the crass predict gives an ontward normal to the body surface. The 
next step is te define the plane of the element by determining the 
averages of the coordinates of the criginai four corner points. These 
points are then projected parallel te the normal vector into the piane of 
the clement to give the corners of the plane quadrilateral. The corner 
points cf the quadrilateral are equidistant from the four points used te 
form the element. Additignal parameters requirec for subsequent force 
calculations, quadrilatera’ area and centroid, may now he calculated, 


The spacing and orientation of the elements is varied in such » way that 
they describe the vehicie shape accurately. Since four puints are used 

to define the plane quadrilateral, the eu,,es of adjacent elements are not 
coincident. This is not important, siace the pressure is caleuiated only 
at the quadrilaterazi centroid. This pressure is then assumed ta be con- 
staat over the surface of the element. 


The plane-quadzilateral surface description method is not as elaborate 
28 some of the other methods. It is important, however, to note thas 
the simplicity cf the method permits the use of conventional crosa- 
sectional drawings in data preparation ine surface slopes required) and 
the use of semiautomatic data-reading techniques. Also, as hae been 
illustrated in Yolume I, computer-generated pictures are uged in check- 
ing the geometric data for errors. 


Parametric Cubic 





A second technique for describing three-dirmensional curved surfaces is 
also provided within the program. This is a mathematical surfece-fit 
technique and is identified as the Parametric Cubic Method because of 
the general type of equatiors used, 


Several different mathematical surface-fit techniques are deccribed in 
the literature. The one used in this pregram was adopted from the 
formulation given by Coons of MiT (Reference 20), In this methed a 
vehicle shape is also divided into a number of secticns or patches. 
The size and location of each patch depends upon the shape of the 
surface. 


Tne basic feature of this method is that only the surface conditions at the 

patch corner points are required to cornpletely describe the surface en~ 

closed by the beundary curves of the patch. The basic problem, bow- 

ever, is the determination of ail the infermation required at these corner 

points, i.e., the surface equation requires corner point surface deriva~ # > 
tives with respect to the parametric variables rather than the X, Y, Z 
coordinates, This has been solved by the use of additional points along 
the boundary curves as will be described Jater. 


i6 


In the following discussions we will use the geometrical representstion 
of a surfane patch as illustrated in the figure below. ; 
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Since the basic surface-fit equations and their derivatives are presented 
in Reference 20, they need be only reviewed briefly in this report. 


The X, ¥, Z coordinates of a point on the surface are related to the two 
parametric variables u and w. Thus, a surface in space is mapped 
into the u, w unit square, The basic problem is to find the position 

{(X, Y, Z) of a point (u, w) in the interior of the section surface, The 
general procedure is to first find relationships for the four boundary 
curves, These are defined as third-order polynomiais in terms of the 
parametric variables. The points on the boundary curves correspond- 
ing to u and w (0, w and u, 0, etc.) are then calculated, A general 
surface equation is used to caiculate the pronerties at the point u, w. 
This equation uses blending or weighting functions to properly introduce 
the influence of each of the related boundary-curve points and the four 
corner points. The blending functions also ensure the continuity of the 
slgpes across the boundaries hetween adjacent sections. 


There are several methods for calculating the direction cosines of the 
tangent vectors required in the calculation of the carner~point deriva- 
tives, Most require the specification of additional surface-boundary 
points, some of which may lie on the extensions of the boundary curves. 
The derivatives must be calculated, since it would not be practical to 
measure them directly from drawings. The method in this program in- 
volves the use of circular arcs through it.ree boundary-curve points, the 
middle cne being a corner point. 


The first step in the cemputational procedure is to determine the equa~ 


tions for the cubic boundary curves. The equation used is given by the 
following relationship for u >= 45. 
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#0, w) = Aw? + Bw" + Gw+D a 
where 
6X, OX, 
A= 2 [x, (0.0) - x, (0,1)] + =-4 (0, 0) +8 40, 2) 
Be 3 i, {G, }} - x, (0. 0)| ~ 255 (6, Qj - “aq (8, 1) 
5x, ‘ 
Gs a (0, 0) 
Ds x; (9, 0} 
Similar equations are needed for the other three boundary curves with u = 1, 


w = 0, and w = tI. 


The miseing items required for the solution of the above equations are the 
derivatives 


° 


8%. ax 
He (0,01, (0, 1), ete 


in the Arbitrary-Body Pragram these are determined by passing a cir- 
cular arc through three points, the middie point being the corner point 

itself, For completeness, the development of this method is presented 
and the sketch below is useful in following the derivation. 
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This sketch is a view of the plane of the circle with U,, as the base coordinate. 
The vectors Tyr T3s and T, are tangents to the curve “at the points 1, 2, and 3. 


The tangents make che angles 5,, 65, and 6. with respect to ee The chord 
lengths make the angles Ey and€, with respect to the vector U3: 


One of the properties of circular arcs is that the chord angle is the average 
of the two tangent angles. 


£2 re BE: Pes ite 
1 - 2 2 2 3 2 


For the coordinate base selected (U,.}, ¢, = 0, therefore, 
13? 2 


a = ~d, and Gs = €) + €5 


The tangent vector at point 2 is then given by 


qT = ces 85 Uis + sin 5,0) ) 

— L = 

U ae 2 » Ly3 is chord vector between points | and 3, 
is [By 


To determine Uw the binomial Uan must first be found 


Yan = big % hip 
U 

Un 2 pol (unit vector) 
Pan 

Un = Yan * Yi3 


The radius vectors (X, Y, Z) fer the three points are 


r, = Xit Yii+ Zk 
tr, = Xi + YJ 4 Zk 
by = Xiit ¥sj+ Z,k 
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Lin, = Fg ° Ey (My ° XPit (¥,- Yo) i+ (2, - 2) 


Ly, = Fg- 8, = (Ky- X14 (¥,- Yy) 9+ (Z,- 2))k 
and the chord angles 
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Fo convenience we will use the shortened notation: 
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; etc. 
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Ue = &,it+m,jtn,k 
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ny (my, n,~ m™m, n) - aT (2) m, ~ i, m, )| j 


+ [m, (m, va 7m, n,) + 4, (2, n, - &, n,) hi 
U, = &,f+ my j+a,k 


And finally we obtain the tangent vector 


R; = (£, cosé + 2, sin 6) i+ (m,; cosé + my sin 6) j 
= (n, cosé t+ ny sin 8) k 
where 
1 i 1 Lis 1 In 
L -[(x,-x)240r,-¥)2+@,- 27] 
13 (late : aera 3 "1 
4G “[, ar lr Mr ely a la m,)| 


my = ~[ ny (my ny - m, nj) + 4, (4; m2 ~ 2, m,)| 
My = [m, (m, n,~ m, n,) + on oat n, - a, a) | 


and 


‘ 
12 ss 12 


2 2 1/2 
po - ~ / 
L = (X5 Xx) + (Y5 Y,) + {Z | 


2 
12 - 2)) 


2 
The final end point derivatives are then found from 


Ue ae OS: 
jw -” i bw 
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where oe 
3s : ‘ Xs 
na the boundary length since Aw = 1 an the unit square patch 
t= NB - 1 i : : 1/2 
aSre es ae eggs ty pga ep 
T= 2 
I = 2 at the starting corner point 


ron) 
t 


NB - | at the final point on the boundary curve 


NB = number of points input on the boundary with one point 
extending off each end of the boundary curve. 


Once “he boundary curves are found the values required for the general 
surface equation can be calculated. This equation is given below, 


X tu, w) X,(9, w)F (u) + X(1, w)F (a) + X,(u, 0)F (w) 


+ 


X tu, 1)F , (w) - %,(0, OF (u)F tw) 


X,(0, 1)F (uF (w) - X,(1, OF (U)F (w) 


X.(1, 1)F, (u)F, (w) 
where the terms Bo and Fi are blending functions given by 


Fy (u) 


3u" - 2u° Fy (u) = Sw" - oy? 


Fv 


1 - Fy (u) F ow) ae F (w) 


The program does not use the parametric cubic geometry data directly in 
the pressure calculations. Instead, the parametric cubic data are used in 


creating surface elements by a systematic variation of the parametric 
variables w, and u. 


One advantag* of the mathematical surface-fit technique over the plane- 
di*tributed-element method is the smaller number of surface points re- 
quired to describe a shape. However, additional points are required on 
the boundaries to determine the required corner derivatives. This method 
is not as adaptable to semiautomatic data-reading techniques, since the 
organization of the required input data is more complex. The accuracy 

of this method depends upon the distribution and crientation of the surface 
sections, just as the plane-distributed-element method depends upon the 
distribution cf the elements. 
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Auxiliary Geometry Methods 


The program has two other geometry generation options. These are (1) the 
ellipse surface generation, and (2) the slab delta geometry generation. In 
each of these methods the program caiculates the necessary Y-Z coordinate 
data at specified X stations. In the ellipse generation simple ellipse equa- 
tions are used; in the slab delta geometry generation, similar equations are 
used along with the necessary equations for the bottom and top parts inboard 
of the nose and leading edge. 


Control Surface Geometry 
The geometry data for a control surface flap are input to the program in the 
undeflected position. The methods used in transforming these data to the re~ 


quired deflected position are outlined in the following discussion. 


The coordinate system used in these derivations is shown in the drawing 
below. 





The general procedure involves a coordinate shift and an appropriate rota- 
tion to a hinge-line centered coordinate system such that the new Y-axis 
(Y,,) lies along the hinge line. For » and $ equal to zero and with the flap 
surface normal in the negative z-direction, the hinge-line centere’ coordi- 
nate system has the same directions as the body-axis system. The corner 
points, centroid, and normal vector (direction cosines) for each element 
of the flap are transformed into this system. Since the flap is a rigid body 
this information is independent of flap deflection and the hinge moment 
factor (moment per unit normal force) need only be determined once. How- 
ever, the force magnitude is a function of the deflection anyle and requires 
having the geometry of the deflected flap in the vehicle-centered coordinates. 
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The coordinate system shift is given by = 
ase ie 
= X- “HL, 
Yio=s ¥-Y¥ 
HL, 
25 Z-2Z 
HL, 
Where 


{ ) is to point 4 on the hinge line 
HL, 


= 


The new coordinates of the flap in the shifted and transformed coordinate 
aystern are given by 


‘ 7 
X5 a 
oe eae y! 
2 Z! 
Where 
ele [ol[+] 
le] - [9] 
cosy sinw 0 
[>| = |-sin) cosw 0 
0 0 1 


wf = rotation about the Z' ~axis 


1 0 0 

@ = aA ° 
[ | 0 cas? sing 
0 -sind cosd 


$ = rotation about the X',-axis 
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The final rotation to the deflected position (6 is the control surface 
deflection) is given by 


! t 
t § 
Mig “2 lo ST Ee 
‘ 
205, Zo Zz 
where 
cosé, 0 - sind, 
[ae | = 0 1 0 
sin 6, 0 cos 6, 


The coordinates of the deflected flap are then transformed back to 
vehicle centered coordinate system, first through the inverse rotation 


t ~ 


t Ord 
& ‘ eer eee 
i 
Y de = [= | 080 
z,' Ze 


and then by the coordinate shift 


' 

Xs, = Xe, + Xu, 
' 

Y6o = Y 66 + Yuta 
' 

Z be = Zé + ZHL4 


The rotation angles are defined for a right-handed system and are 
found from the relationships 


.1 ( *Hu, > *HL -1 { 4HU, ~ 4HL 
wos sg” ( “Ha ~ “is and $= = sin7! page ee 
XY 
where r 2 Ip 
Lyy = ce - Xyn,4) + (YHb,) 
and ; 3 Y> 
YZ. = | xy + (2Hr, - 2HL4) | 
a 


A check is made in the program and if Yur, < YHLg then the yaw 


wv = @ - w to position the hinge line in the 


rotation angle is set to 
proper quadrant. 


The third rotation angle & is, of course, specified for a given prob-~ 
lem. Jt should be noted in the present approach, that the coordinate 
system is rotated through the angle 6,, positive in the right-handed 
sense for the system defined. Relative to the physical problem, posi- 
tive 6, corresponds to a flap deflection into the flow. 


The hinge moment factor (HMFCT) is simply a function of the elernent 
geometry and location, and is defined as follows. The total moment 
of an element ie (considering only inviscid forces) 


wa 


My 


tt 


- (Ry x F} = P(R'g x Np) AREA 
where 
is the radius vector to the elernent centroid, 


is the net surface pressure, 


and AREA is the element area. 


The hinge line moment is just the ¥'g- component of the total moment; 
Myr = My = j'9 ®° M'g = P (HMFCT) 


where 
HMFCT 


(Z'g Ny'g - Xo Nz ()) AREA 


Once the deflected flap is properly oriented in the vehicle centered co- 
ordinates, the force on each element and hinge moment are determined. 


Graphics -~- Picture Drawing Program 


The perspective drawing computer program is an important component of 
the Arbitrary-Body Hypersonic Force Anaiysis System. Its use in this 
system is in providing graphical perspective drawings of the geometric 
description input to the arbitrary body force program. * The purpose of 
these drawings is to allow the engineer to detect errors in the geometric 
input data to the arbitrary body force program. 


As explained previously, the body shape is defined by input sets of points 
in turee-dimensional space. A grouping of four surface points is used to 
describe a surface element. An organization of a large number of related 
elements forms a body section and a number of sections inay be used to 
give a complete description of the shape. The equations required to pro- 
duce the perspective pictures are derived in the following paragraphs. 


*The drawings made by this program are noi true perspective drawings 


but are a limiting case where the vanishing points have been moved to 
infinity. 
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Each point on the surface is described by its coordinates in the body 
reference coordinate system. 


xX 
Y 


The body reference coordinate system is assurned to be a conventional 
right-handed Cartesian system as illustrated below. 





x 


To create the perspective drawings illustrated in this report each 
surface point on the body must be rotated to the desired viewing 
angle and then transformed into a coordinate system in the plane 

of the paper. With zero rotation angles the body coordinate system 
is coincident with the fixed system in the plane of the paper. 
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The rotations of the body and its coordinate system to give a desired 
viewing angle are specified by a yaw-pitch-roli sequence /, 6, 0). 
This rotation is given by the following relationship: 
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Since each point on the surface it given by its coordinates in the 
X, Y, Z system,its position in the fixed ccordinate system (Xo 
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zZ.) may be found by reversing the above process, 
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If we carry out this operation we obtain 


cos8cosl -sin? coedétsind« osh sind 


=jcos6sinu cosh) cosétsinOsins sind 


~sing 
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x, =X{cosécesy) + Y¥(-sinycas$tsinscos}sing} + ZlsingsindtsinGcosbe0sd} 
ae =X{cosGsinu} + Yicospcosétsingsindsind} +Z{-cosdsinptainésintcosd) 
Zz, =X(-sin6} + Yicos8sin $ + Zices8cozd} 


We may now use these last two equations to transform a given point on 
the body {X, Y¥, Z}) with a specified set of rotation angles (¥, >, 9) inte 
the plane of the paper (the Yo Z system}. With the SC-4020 library 


subroutines it now becoma2s a simple matter to plot these data and to 
connect the related points with straight lines. 


In the surface fit technique used in this program and described in 
Reference 22, each input element is replaced by a plane quadriiateral 
surface element whose characteristics.are used fer ali subsequent 
calculations. These characteristics include the area, centroid, and 
the direction cosines of the rurface unit normal. The surface unit 
normais may be transformed through the rcquired rotation angles 
just as was done for the individual points. The resulting value of 

the component of the unit normal iu the x, direction {out of the plane 


of the paper} may be found fram the following equation. 


ee nleosGcoosin, (-sinucosbfsinGeostamd)ta, (sindetnetainecospcos¢) 


where ane ny n, ere the components of the surface init normal in the 


vehicle reference system. 


if nh, 18 positive then the surface element is facing the viewer. If ny, 

° ° 
is negative the element faces away from the plane of the paper. This 
result is used in the program to provide the capability of deleting most 
of those elements on a vehicle that norrnally could not be seen by a 
viewer. The resulting picture is thus made more realistic and confusing 
elements which are on the back side of the vehicle do not appear. No 
criterion is provided, however, for the deletion of those elements that 
face the viewer but are blocked by other body components. This may 
be accomplished by a proper selection of viewing angle or by a physical 
deletion of the offending section from the input data. 


COMPUTATION OF VEHICLE FORCES 


Calculation of Lecal Flow Conditions 





In the preceding derivations we have converted the input element into 
a plane quadrilateral element. The quadrilateral is described by its 
area, the coordinates of the centroid of the elerment, and by the direc- 
tion cosines of the surface unit normal. In the force calculation 
methods we mast also know the angle that the element makes with the 
free-stream velocity vector (the impact angle). This angle changes as 
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the vehicle attitude {angle of attack and yaw angle} changes. The 
imspact angle may be found from the following relationship: 


6 = Fle «6 
eV 
COS Qo eee 
la] {71 
where 
R is the unit normal cutward from the surface with 
direct on cosines ny , Ry, Nz 
Vv is the local velocity vector with directien cosines in the 
vehicle coordinate systern given by V,, Vy Ve 


The direction cosines of the unit surface normal are given by the quadri- 
lateral calculations. The value of the local velocity vector V depends: 
upon the vehicle attituge with respect te the free-stream direction and 
its angular rotation rates, and is derived in the discussion below. The 
rotation directions are consistent with the conventional stability body- 
axis system. The coordinate systern, however, is changed to be 
consistent with the geometric descviption system discussed previously. 





wnere 

rolling velocity 
pitching velocity 
yawing velocity 

total angular velocity 


roe 2 Out) 


Rud 


The movement of a given element of the vehicle with respect to the free~ 
streara depends upon the vehicle rotation rate and the position of the 
element relative to the rotation center. The radius vector from an 


arbitrary reference point on the vehicle to a point on the surface is 
giver. by 


Fos (x~x)irly-y) it (e-4)k 
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where ¥), Yg» %, 18 the moment reference point (center of gravity). 


The total angular velocity is given by 
$2 PT.6F-Re 
The free-stream velocity vector is given by 


wo food 


¥= Gait Vay J Veo, 


The free-stream velocity components are given by the following 
relationships for a conventional yawepitch sequence. 


- V,, Cos @ cosR 


Ox 
Vooy, = V,,31n B 
Von = Yy cosp sing 
where 
ae = angle of attack 
8 = sideslip angle 


The total velocity vector relative to the surface element is obtained 
by combining the above relationships as follows: 


Vv 7 Vo-axr 
a 


The local velocity vector therefore becomes 


Vo |¥a + [a ten) Rory) 3 
x ae 
+ 1 + IR (x-x,) + P (z-z j oti 
4 
+ va > Piry) toe) 
{ Zz 
or Vos V_itvV j+v_k 
x y Zz 
where 
v= = [Q (z-2,) - R (y-y,)] 
vy = vVo+ [R (x-x,) + P (z-z,)] ‘ 
at F [P (y-y,) + Q (x-x,)] 


Zz 


The total local velocity is given by 


oe ey ere eee 
Viecal = Vise + Vy + Ye 
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The coaventional surface impact angle is then given by 


6 = 0/2 - coat (Ee 727) 
Viocal 


where 


Ny s My» 2, are the cutward surface unit normal direction 
. cosines 


Vehicle Coefficients and Derivatives 


In the program force calculations the pressure on each element is calculated 


" completely independent of all other elements (except the shock-expansion 


method). If the vehicle is rotating the local pressure coefficient must be 


corrected back to free stream conditions. This is accomplished by the 
following relationship, 


2 
sr fe Migesi 
GC = G <a 
P Plocal od 


in the arbitrary kody force calculation program the pressure coefficient 
on each element is calculated completely independent of all other elements. 
The contributions of all the elements are then summed to give the total 


ferce components. The basic relationships to accomplish this are given 
below: 








axial ferce Cc = A Cc n AA 
A 5 p x 
ref 
l 
side force Cc = —- xc n AA 
¥ Saf Poe 
normal force Cy = -= ye n AA 
ref P 
where 
AA = element area 


Note: The minus sign is needed because of the sign convention 
on Zin the body coordinate system, directed upward. 


The moment coefficients are obtained by a summation of the 
component forces multiplied by the distance from the element 
centroid to the reference moment center. 


32 


ep 


~~ Fe 


vv 
rolling moment CG, = 2 Ey a + S (Gry ra) 


‘ : a: F x Zz 
pitching moment Cc = S(Cy =) + Sic m =) 


awi = S x - mA 
yawing moment Cc. ice 5) Sc oe 4) 
where 
b = reference span (lateral and directional 
moment coefficient reference length) 
¢ = mean aerodynamic chord {for longitudinal 
moment reference) 
MVPS > SS distances froin the center of gravity 


x 7 7 xX ate. 
centroid cg’ 


The conversion of the axial force and normal force coefficients to lift and 
drag roefficients requires the following rotation matrices. 


Cr -1 Cy 
! = EK 
Cy [=| rey 
LoL Cy 
where 
1 cosa cosg -sing sing Coss 
[z}° = | cosa@ sing cose sina sing 
-sing 0 cosa 
a = angle of attack (+ nose up) 
8 = sideslip angle (+ nose left) 
Ch = Ca cosa cosep - Cy sing + Cn sine cos@ 
Cy = Ca cosa sing + Cy cosp + Cn sing sing 
Cy = -C, sine + Cn cosa 
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The vehicle static stability derivatives in angle of attack and sideslip 7 
are calculated by the method of small perturbations as indicated Se 
below: 
c = (Cy), tad (cq) 
Ao Aa 
Cc ~ {Cc 
Cc = ( N \ tAa ( N)o 
Ny ha 
cS Cc 
oS = ( mo tae ( m), 
My 
Cc CG 
Cc = ( Y (Cy), + AB ~ ( vp 
¥6 AB 
The damping derivatives due to vehicle rotation rate are given ina 
similar manner 
r 
c - Cin qtAg ~ (mm) a | 3 
m, AG 2V 
etc, 
The control surface derivatives are also calculated by the methced of smali 
perturbations. 
Cc - (Cc 
ae ( a +A6 ( L}, 
Leg “ As& 
Cc ~{C 
ae ( m), + A5 ( n), 
mM ¢ PAN) 
C ~(Cy 
oe ) 6 +A6 ( } 
ie Ad 
Cc -{C 
: ( v), as y}, 
Cy: 7 As f 
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Inviscid Force Calculation Methods 


Many of the pressure calculation methods used in the analysis of high-speed , 
shapes are listed in Figure2 . An attempt has been made in the prepara- 

tion of this figure to indicate the interrelationships of the methods (the in- 

formation can, of course, be organized in many different ways), Some of 

these methods are more applicable to the arbitrary-body problem than 

others, 


The method of characteristics is the eventual ideal approach for the cal- 
culation of forces on three-dimensional shapes at high speeds. It will re- 
quire starting solutions for three-dimensional blunt bodies of acbitrary 
shape. The development of a method of caiculating three-dimensional 
boundary layers would permit the use of an iterative process to account 
for the viscous~inviscid interaction, Although this approach has been used 
for some very simple shapes, the complete solution for arbitrary shapes 
is some time away. Significant progress is being made in the solution of 
the inviscid flow fie id by the method of characteristics. However, present 
mathematical techniques and digital-computer size and speed capability 
prevent application to typical preliminary-design problems. Application 
must be reserved for simple shapes or important detail design applica- 
tions where very large computer times might be acceptable. 


Many of the other methods shown in Figure 2 would be useful force- 
calculation methods for inclusion in an arbitrary-body system. The selec~ 
tion of the proper method in a given application depends upon the vehicle~ 
component shape and flight condition and must be selected by the engineer 
on the basis of his knowledge and experience in the use of each method, 
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Figure 2. Pressure Calculation Methods 
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The arbitrary body force computer progra.n contains a nurnber of optional 
methods for calculating the pressure coefficient. In each s.ethod the only 
geometric parameter required is the element impact angle, 6, or the change 
in the angle of an element from a previous point. 


Before the program calculates the pressure on each surface element, it 
checks to see if the element is facing the flow (in an impact region) or 
facing away from the flow (in a shadow region}. The methods to be 
used in calculating the pressure in impact and shadow regions may be 
specified independently. A summary of the program pressure options is 
presented below. 


PRESSURE CALCULATION METHODS — MARK ill MOD 0 PROGRAM 


Impact Flow Shadow Flow 

1. Modified Newtonian 1, Newtonian (Cp = 0} 

2. Modified Newtonian+Prandtl-Meyer 2. Modified Newtonian+Prandtl-Meyer 
3. Tangent wedge 3. Prandtl-Meyer from free-stream 
4, Tangent-wedge empirical 4, OSU blunt body empirical 

5, Tangent-cone empirical 5, Van Dyke Unified 

6, OSU biunt body empirical 6, High Mach base pressure 

7, Van Dyke Unified 7, Shock-expansion 

8. Blunt-body shear force 8, Input pressure coefficient 

9. Shock-expansion 9. Free molecular flow 
10. Free molecular flow 10. Boundary-layer induced pressures 


11, Input pressure coefficient 

12. Hankey flat-surface empirical 

13, Nelta wing empirical 

14, Dahlem-Buck empirical 

15, Blast wave 

16. Modified tangent-cone 

17, Boundary layer induced pressures 


Since most of these methods are adequately discussed in the literature they 
will be reviewed only briefly in this document. The blunt-~body shear force 
and the boundar,~-layer induced pressure methods are discussed in detail in 
the section describing Viscous Force Methods. 


Modified Newtonian 


This method is probably the most widely used of all the hypersonic 
force analysis techniques. The major reason for this is its simplicity. 
Like all the force calculation methods, however, its validity in any 
particular application depends upon the flight condition and the shape 
of the vehicle or component being considered. Its most general ap- 
plication is for blunt shapes at high hypersonic speed. The usual 
form of the modified Newtonian pressure coefficient is 


C =  Ksin“s 
P 
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In true Newtonian flow {M=o, Y= ]} the parameter K is taken as 2. 
In the varioug forms of modified Newtonian theory, K is given values 
other than 2 depending on the type of modified Newtonian theory used. 
K is frequently taken as being equal to the stagnation pressure co- 
efficient. In other forms it is determined by the following relation- 
Ship (Reference 23}. 


Cc 
K _ Prose 
7 2 
sin 6 
nose 
where 
Cc = the exact value of the pressure 

Prose coefficient at the nose or leading 
edge 

Boas = impact angle at the nose or leading 


edge 


In other work K is determined purely on an empirical basis. 


K = fn (M, @, shape) 


When modified Newtonian theory is used, the pressure coefficient in 
shadow regicus (5 is negative) is usually set equal to zero. 


Modified Newtonian Plus Prandtl-Meyer 


This method, described as the blunt body Newtonian + Prandtl-Meyer 
technique, is based onthe analysis presented by Kaufman in Reference 
24. The flow model used in this method assumes a blunt body with a 
detached shock, followed by an expansion a~ound the body to supersunic 
conditions. This mathod uses a combination of modified Newtonian and 
Prandtl-Meyer expansion theory. Modified Newtonian theory is used 
along the body until a point is reached where both the pressure and the 
pressure gradients match those that would be calculated by a continuing 
Prandtl-Meyer expansion. 


The calculation procedure derived for determining the pressure co- 
efficient using the blunt body Newtenian + Prandtl-Meyer technique 
is outlined below. 





1, Calculate free-stream static to stagnation pressure ratio 
ae at 
p 5 ae 2YM,,” = (vy - 1) Y-1 
P = St —_—_— a ee aR a! 
Po (Yt 1) MQ" Ved 


ies) 


Assum? a starting value of the matching Mach number, M 
(for Y= 1.4 assume M, = 1. 35) q 


3:7 


ese, wer ee 


10, 


Calculate matching point to free-stream static pressure ratio 


pases ep 
Y-i 


Assume a new matching point Mach number (1.75) and repeat 
the above steps to obtain a second set of data. 


With the above two tries use a linear interpolation equation to 
estimate a new matching point Mach number. This process is 
repeated until the solution converges. 


Calculate the surface slope at the matching point 


se .~ Q-P 
sin by = i. Pp 





Use the Prandtl-Meyer expansion equations to find the Mach 
number on the surface element, M; 


Calculate the surface pressure ratio 


apes 
ies. x fs pba | ne 
Po an 
where 
uP is provided as an empirical correction factor 
Ps is the pressure on the element of interest 


Caiculate the surface to free-stream pressure rafio 


P \ 
Hot os) { 
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il. Calculate the surface pressure coefficient 


Pp 
C = ees So ee | 
Fe Via Vrs 


The results of typical calculations using the above procedure are 
shown in Figure 3. Note that the calculations give a positive pres- 
sure coefficient at a zero impact angle. As pointed out in several 
references these results correlate well with test data for blunt 
shapes. However, if the surface curvature changes gradually to 
zero Slope some distance from the blunt stagnation point the pres- 
sure calculated by this method will be too high. This is caused by 
characteristics near the nose intersecting the curved shock system 
and being reflected back onto the body. If the zero slope is reached 
near the nose (such as in a hemisphere or a cylinder) this effect has 
not had time to occur. 


Tangent-Wedge 


The tangent-wedge and tangent-cone theories are frequently used to 
calculate the pressures on two-dimensional] bodies and bodies of 
revolution, respectively. These methods are really empirical in 
nature since they have no firm theoretical basis. They are suggested, 
however, by the results of more exact theories that show that the 
pressure on a surface in impact fiow is primarily a function of the 
local impact angle. In this program the tangent-wedge pressures are 
calculated using the oblique shock relationships of NACA TR-1135 
(Reference 25). The basic equation used is the cubic given by 
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Figure 3. Biunt Body Newtonian + Prandtl-Meyer ; 
Pressure Results 
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The roots of the above cubic equation may bé ebtained by u Sing the 
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= g0jution procedure (see Reference 26} 


indi dicated 


2 - p/3” cos {wf/2} - b/3 
-2 J= p/? cos {w/3 + 60°) ~ b/3 


-2 J- p/3 cos (u/3 ~ 60°) - b/3 


R, = y, 7 b/3 
R, = Y> ~ b/3 
R, = ¥3 7 b/3 


roots of the reduced cubic equation 


pe 
ag ee 


21b/3)° = a: 


a ee 


2 ¥~-{p/3) 


4 : : 
sir 8, = roots of the cubic equation 


The smallest of the three roots corresponds to a decrease in entropy 
and is disregarded. The largest root is also disregarded since it 
never appears in physical actuality. 
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‘For small deflections, the cubic golution becomes very sensitive to 


numerical accuracy; that is, to the nurnber of significant digits car- 
ried. Since this is dependent on the particular machine emplo ;ed, 
an aiternate procedure is used. 


Whan the flow deflection angle is equal to or less than 2.0 degrees, 
the following equation is used instead of the above cubic relationships 
(Reference 27}: 





3in 


2 ] . ¥t) 6 
. ‘ 2 
M 


Once the shock angle is obtained the remaining flow properties may 
be found frem the relationships of Reference 25. 





‘GMS sin” 8. 
density = Po = p oo 
F M” sin 8, +5 
7(M° sin” 9, - 1) (M° sin? 9, + 5) 
temperature = 15: 7} 3 
[ 36 M~ sin 6, 
7™M® sin” 9, - 
pressure ia See 
coefficient = Cy = eae ern 
0.7 M 
where 
{ 32 = conditions behind the shock 


Oblique shock detachment conditions are reached when 110 solution 
may be found te the above cubic relationships. Under these conditions 
the program uses the Newtonian + Prandtl-Meyer method for continued 
calculations. 
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Tangent-Wedge, Tangent-Cone, and Delta Wing 
Newtonian Empirical Methoa 


The tangent-cone and the tangent-wedge Newtonian empirical methods 
used in this program are based on the smpirical relationships derived 
below. 





For wedge flow 





sin by 
sin 6_= ———— ee 
s - - 
(1 - €) cos (6. 5.) 
where 
p Y=) 2 
€ — —_— = I +- ooo: 
p Y+ 
2 d ive ipa? 
ns 
For cone flow (thin shock layer assumption) 
sin 5. 
Sin @, = 7 perry ane 5 ; 
~ 3) cos (8. - 6.) 
In the lirnit as M--x, € = & x teed and cos (@_ - 6) = 1 
; lim Yt s 
Therefore 
wedge cone 
Yt : ! 
oinQ - Seta i ee, sin @ = AYE) sin 5. 
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These limiting expressions fer @ may now Le compared with the 
data of TR-1135 (Reference 25) at y= 7/5 using the following 

similarity parameters. The exact equations contain three vari- 
ables — Gg, 6, and e«. Noting that for Y= constant, ¢ = fn(My5) 


only, the preceding equations may be rewritten in the following 
form: 


wedge 


cone 
M sin 8. M sin 6 
M = OS Moos 
ns {l - €} cos (8. - é) ns 


€ 
fl - 3) cos (6. - 6.) 


The parameter (£ - 6) is approximately constant and independent of 
M except near the shock detachment condition. The equations es-~ 
sentially contain only two variables, M__ and M sin 6. These are 
used as coordinates to plot the data for Wedge flow shown in Figure 
4, A similar plot could be obtained for cone flow. From the figure 


it is seen that the data are nearly normalized with the use of these 
coordinates. 


For rapid calculations we need relationships for M 
M sin 6 that satisfy the following requirernents: 


5 28 4 function of 
li The effect of shock detachment is neglected 
2. At M siné = 0, M., 7} 


3. The solution asymptotically approaches the M =2 line 


dM 
ns : a 
4, Have the correct slope, aM sing ee M siné = 0 


These conditions lead to equations of the following form 


Ry 
-=~ M' 
wedge M._ = K Mite ~ 
ns Ww 
. Yel 
Ky * 2 
“K, M! 
= t 
cone Mos Ko M'+e 
where 
M! = M sin § 
K. 2 24¥ 4 Y/(y + 3) 
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Figure 4, Wedge Flow Shock Angle 
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These expressions are compared with the data of TR-1235 in Figures 


5 and 6, The cone data are also shown in Figure 7 with the same 
scales asin Figure 4. 


The pressure coefficient may now be obtained by the following 
relationships for a wedge and cone respectively. 


c 


$ 


> = (eq) (Mas* Df Me 


‘?) 
" 


(y= 1) Mys2 +2774 
4(yt]l) My. 


“ei 


[ 
2 sin? ‘ = 


Experimental results have shown the pressure on the centerline of a delta 
‘wing to be in agreement with two-dimensional theory 2t small values of 

the similarity parameter (M'< 3.0) and with conical tlow thesry at higher 

values, The previous expressions derived for wedge and cone flows have 


been combined to give these features, The resulting relationships are 
given below. 


K 
' “Ke . =) M! 
Mas = Keo Mi + e 
For Y = 7/5 
Mas = 1.09 Msin§ + e 70-49 M sin6 


The similarity parameter relationship for pressure is 


Mm’ ©, : (34) (m,.° 7 1) 


The shock angle and pressure coefficient calculated from the above 
equations are compared with the experimental results (Reference 28) 
in Figuies & and 3, respectively. 
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Wedge Flow Shock Angle Empirical Correlation 
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Figure 7. Conical Flow Shock Angle Empirical Correlation 
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o 


ws 


GN te 


h 


i 


vation 


Ate Ree ae ee. et encs 


5 fas r 
4 ape SUE OI RAE het ak CT ORR OR UA Rg ett mehonanaity kamen CR WA OR me h 





ips: ake at Oa Sat aah Se de 





EEE Ee ECT EES Ee 
HEA AS aah rH ae EHH ry seageee mae 32 me se SSHMMRIN SLES RERS: Per 
mils Eee faite disse Sailisias pratt Pai cece a PETER 
BREECH cna HH HIRE EEE IR un EER 
le Pee R ee eet deaaliai uiceesees He eeSsuuii HH sipieseasee 
HEEL csreetigliad Fala eteeaam gue HESsEe 
SHEE EE | ade ilttcas sites es Eee EH sian: iisseeugos 
sciuantts ERS ESRE MTHISAT SS SERS SISSGA SSE ES HTIHTSLHESOSNTIHESESDoSEOS 

: Hale REE Hat emg eee ee ee ae Hstsaseeee 





¢ 
rs 
aubea 


LL 
HE 
HH 
rau 
eG 
yi 
Seer 
ekantnet py 
an 












7 Sarat) eS anTicesee 
| a Site ciessgetae aa EEE Ee 
atte erect ae a fsa 
HES | dine HE iti ti Hee TEEPE 

tI aie ia Sth hast tI nan Ht THT ee tT TP | 
Hees =f srse — Siassd satis Wasi ease ceeag Oe aE STORES Ha TTTETEE SOSA TT RTTTaE ES 

TP ©, + 4 Ho HH poste ae oa. 

iitsta te HERE 





ene 9 eben whee 8 aa eae 


Srreneaat faa SS 
a : a 
cee 


a 


: : =o 


Fea 
| 


rag pegged 


100.0 Ferra 


HEE HN ett = os iatatenwan=sae cl (usenienl 
EEE ee 
S 


a 
SB Seeesenees yuiiisiscees! : 
cen Hee iutgreeaae = 


Euea ga Ginter eceteeeSEe 
geelevte itt tapreseae 





rt 2 ay a OP 


nage elk 
ane % a ices a ae pirisaeece 
Useese: sifcse Pee! SHES STN HI eRe 


L 
i 


ti i 


srt 


oe samerer wenen ae me oar 


ae oso 


Gatti iices 


| a Opes mae eae 
Pen EE 


S 


ar 


o 


MSing ~ Mt 


Delta Wing Centerline Pressure Coefficient Correlation 


Figure 9, 


51 





ERAT EERO ae 2radvnchntteereetesonte Fc rvetunentsent aden ay {ramet pemRE HL NR AMETREENEL ACTIN AN YATE UFO STTSCAANIE ES SY ANE DTI TT RY PS SAI ON AP it et NE 





fe aT aye. 
POR, Nee aR a? 
a 


Oe Rae 


C5U Biunt Body Empirical Method 


The OSU (Ohio State University) biunt body empirical equation 

desecrives the pressure distribution about cylinders in supersonic 

flow, The eauation was presented in Reference 29 and was stated 

to match "all the data obtained cn the cylinders in the present test 

series With a maximum deviation of 2.5 percent." The expression used is 








—= 0.32 + 0.455 cas®@ + 0.195 ces Z8 + 9.635 cos 38 - 0.005 cos 46 


my SALA ed carr Miia le Pr Sr, Sakai ne aelt Ses 
. We! . set 
Ct & : 
, 
y | 


8 = perighera! angle on a cylinder 
-(= 0 at the stagnation point) = 90° -& 


ee 


Pp; = surface pressure 





= total pressure rise through normal shock 


ie aL 
w. f 


The pressure coefiiciert is calculated from the relationship 


3 P 
> ~ 4 FBy c\. 
Cc = tl — 4 (Sra?) 
P Pe, | 


x 
4 


Pe 
— = xomeel 
PL z 


i” 
. K = stagnation pressure coefficient = C 
a _ Pstag 
gE 
@ pF freestream pressure 
¥ = ratio of specific heats = 1.4 


Van Dyke Unified Method 


This force calcularicn method is based on the unified supersonic- 
hypersonic smail disturbance theory proposed by Yan Dyke in 
RezZezence 39 as applied to hasic hypersonic similarity results. 

The methed is wseful for thin rrofiie shapes and as the name implies 
extends down to the supersonic speed region. 
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The simiuarity equations that form the basis of this method are 
derived by manipulating tie cbhiique shock relations for hypersonic 
flow. The basic derivations are shown on pages T53 and 754 of 
Reference 33, The result obtained for a cormpression surface under 
the assumption of a smal: deflection angle and large Mach nuinber 
is (hypersonic similarity equation}. 





f #23 
H / re. 
Cc, x= a rcs a aera! ee 


where His the hypersoa‘c Sirallarity parameter given by M&. The 
contribution by Van Dyke in Reference 30 suggesis that this relation- 
ship will also be valid in the realm of supersonic linear theory if 


the hypersonic similarity pararneter Mé& is replaced by the unified 
i 2 4 

3aupersonic-hypersonic pavhmeeenly MU - r)a, This latter param- 

eter is used in the calculations for this force option in the arbitrary 

body program. 


A similar m2thod may also be obtained fox a surface in expansion 
flow with no leeding edge shock such as on the upper side of a> air- 
foii. The resulting equation is 


2y 
z_2 2 en an 
oc = & 4, ja-Liyp “3 
P yHe | 


io. 


e . 2 \ - “ge 
where again H is taken to belv ni - 1J% in the unified theoi , approach. 


Shock-Expanzion Method 


This force calculation method is based on ciassical shock-expansion 
theory {see Reference 31). in this method the surface elements are 
handled in a “strip-theory" manner. The characteristics of the first 
element of each longitudinal strip of elements may be calculated by 
oblique shock theory, by conical flow theory, or by a Prandtl-Meyer 
expansion, Downstream of this initial element the forces are cal- 
culated by a Prandti-Meyer expansion. 


By a preper selection of the element orientation the method rnay be 
used fur both wing-like shapes and fcr more complex body shapes. 
in this latter case the m2thod operates in a hypersonic shock-expansion 


theory mode. 
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Free Molecular Flow Method 


very high altitudes conventional continuum “ow theories fail and 
one must begin te consider the general macroscopic mass, force, 
and energy transfer problem at the 30ey Surface. This condition 
eccurs when the aiz is sufficieatly rarefied so that the mean free 
path of the moiecuies is much greater than a characteristic bed: 
aimengion. This conditicn is known as free molecular flow and tie 
method cf analysis selected for thie program is described in 
Reference 32. This method was alse used in Reference 19, The 
equations used were taken from these refere;ces and are presented 
below. 


Pressure coefficient 





if2-z ft 
~- } I wo: n b ee sin 4) 
G = ome a § sin § > = 
P g*° | a) ¥ aos iy 
Be De ey te AB | \ 
+ Uf2 ~£ MS sin & +3] tha Ssiné [as extis sin oH 
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Shear force coefficient 





(cos ae [ 4 2 
a7 (S sin 5) + 


Cc = = N\x§ cin & |} + erf(S sin 4) 
om v7rS ' 
} 
where 
S = speed ratio = ¥¥/2 M, 
fn = normal momentum accommodation coefficient (=0.0 for 
Newtonian and = 1.0for completely diffuse refiection) 
é = impact angie 
Ty = body temperature, °K 
T. = free-stream temperature, °K 
x 3 
F 2 “x” 
erf = error function erf (x)= T= e dx 
i G 
£ = tangential momentum accommodation coeificien: 


{= 0, for Newtonian flow and 1.0 for completely 
d’=tuse reflection) 
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The pressare force acts perpendicular to the surface and this uirection is 
readily cbiained since the element normal has already heen determined in 
the geometry subroutines. The shear force acts in the direccion of the 

en 


angential velocity component on the surface and this direction is deter- 
i 


mined by taking successive vector products as foilews. 


The procedure is illustrated in the 
accompaning sketch where the in- 
cide:.t velocity vector is defined as 


Vc Vy it Vy 
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First, a surface tangent vector {T} is defined by the cress product of thr 
normal! and velocity vectors; 
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Then the direction of the shear force (S) is given by the cross product of 
the surface tangent ard normal vectors; 


S = 5 i+S.j+5 7k 


x Y 
where 
s ha bc a Ww 
Sy = Ty Ny T, Ne 
Sy es tT, Ny - Ty Ny 


= J ce N 
37 Ty vy Ty “xX 
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. given by 


“in &sing thé fee molecular flow method the above analysis must he 


‘Phe final components of the shear force in the vehi 


icle axis system are 


a on we steep eE: s “ . oe e 
SHEAR, = (SHEAR) (Sj / STOTAL 


” SHEAR - (SHEAR) (Sy) / STOTAL 
SHEAR, = (SHEAR) (sz) / STOTAL F 
where | = "3 _—s 


SHEAR is the shéar force a3 calculated by the free molecular flow 
equations, — ~_S 


: fas 31/2 : 
; 2 2 z : 
STOTAL 2 {8..°+5.,°+55° - 


carried out over the entire surface of the shape including the base, 
shadow regions, etc. When the free molecular flow method is 
selected, it is used for both impact and shadew region. 


“Enis method of determining the shear direction is also used for the 
continuum viscous. forces discussed in the next section. The plane 
forraed by the velocity vector and the surface normal is referred to 
aa the velocity plane (shaded region in the sketch}, since both the 
‘incident and surface velocity aré in this plane. This definition is 
correct for two-dimensional flow, however, it is only an approxi- 


mation tothe shear direction in the general arbitrary-body case. 


Hankey Flat-Surface Emprical Method 


This. method uses an ernpirical correlation for lower surface pressures on 
blunted flat plates. The method, derived in Reference 14, approximates 


tangent-wedge at low impact angles and approaches Newtonian at high impact 
angles. The pressure coefficient is given by 


C, =i. 95 sin? § + 0,21 cos & sin 


Dahlem-Buck Empirical Method 
This is an impact method that has been derived such that tangent-cone and 
Newtonian results are approximated, respectively, at low and high values 


of the impact angle. The empirical relationships presented in Reference 33. 
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are re 
ee, > in gaye Sh 
for 8 < 22.59 c, = tbs ig (sin 617 
; : {4 cos @ ces 26)7° ~ 
a ee Ae pS 
for 8522.59 | C,= 2.0sin’g 


Blast Wave Pressure Increments — : ; é 


This method usss rcanventional blast-wave parameters to calculate the over- 
pressure due to bluntness effects, Force contributions determined by this 
procedure must be added to the regular inviscid pressure forces {tangent~ 
wedge, taugent-coie, Newtonian, etc.) calculated over the same vchicle 
geometry, The specific blast wave solutions used in the Program were de- 
rived by Lukasiewicz in Reference 34; : . 
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where 
: Cp is the nose drag coefficient. 


d is the nose diameter or thickness 


Xo is a coordinate reference point 


and the coefficients A, B are 


Flow 






Two-Dimensional 0.56 


Axisymmetric 7 0.067} 0.44 
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Modified Tangent: ~Cone Method espe ae  : wi 


This method. originally deecioped for vse on cones with elliptical cross- 


‘Sections, modifies the tangent-cone result by an increment representing | 


the deviation from an average pressure divided by an average Mach number, 
Mote specifically, the following equations are used {after Jacobs, Reference: 
ane ; a in 


‘y 


Ps od, pee Gia. cx Site = Crave : 
awiiere oe C_ _is the surface pressure coefficient ; 


eed 
BS Tes is the conventional sangentr eons pressure coefficient 


is the average pressure coefficient 


= USp¢ AlQSs 4 Ai 5 element area 


; aeoerd a 


— 


Mayg is the average Mach number, defined for an eqaivaletit 
: cone having pressure coefficie: tC, ave’ 


igh Mach Base“Pres ates 


For a-body in high speed: ‘low it might be aisateea that any base regions 


‘would experience tofal vacuum, that is, 
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i? 


Sr 


Cc. s » 
P M 


he 


However, the viscosity of real gases causes sume pressure to be felt in 
base region and experimental data have shown this to be roughly 70% 
vacuum for air. Theyefore, the expression 


has been included in the program. 
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Viscous. Force Calculation Methods 


oe 


The most difficult part in the analysis of an.arbitrary shape is the calénla< 
tion of viscous forces. A detailed knowledge of the local properties and on 
flow history along surface streamlines is required, This combined with. the 
natural complexity of the boundary-layer equations necessitates considerable 
simplification of the problem before solutions can be obtained, An engineer- 
ing approach has been selected that retains the essential characteristics of 
the hypersonic houndary-layer problem. No attempt is made to calculate 

the detailed skin friction distribution on the exact arbitrary shape, but rather, 
the vehicle ig represented by a number of flat surfaces on each of which the 
sheax force is determined, 


The surface streamlines are assumed in the velocity plane and the flow 
history is approximated by the inclusion of.an initial surface, The shear 
force is detetmined fu. both laminar und turbulent flow and may be summed 
over the vehicle for either type, 


Reference temperature and reference enthalpy methods are available for 
both laminar and turbulent flows and, in addition, the Spaiding-Chi method 
with either temperature or énthalpy ratios may be selected for turbulént 
calculations. The surface temperature may be either input or the radiation 
equilibrium value determined. The effect of planform shape, leading edge 
viscous-interaction, and the viscous forces on blunt bodies are also con- 
sidered, 


Skin Friction Geornetry Model 
For the skin friction calculations a geometrically complex vehicle is divided 
up into a number of plane surfaces in a manner which adequately approxi- 
mates the true shape. Leading-edge surfaces and local curvatures are 
omitted, Regions of relatively large curvatures can be represented by using 
@ greater number of plane surfaces, The degree tv which this is done will 
depend upon the complexity of the actual shape and experience of the de- 
Signer, The geometry data for the skin-friction geometry model is prepared 
in the same way as the surface element data used for the inviscid pressure 
calculations and retain their relative location to each other and to the flight 
path. This skin friction modeling technique is best described by viewing, 

for example, a typical high L/D vehicle shown in Figure 10, The upper half 
presents the skin-friction representation of the vehicle which is to be con- 
trasted with the detailed inviscid geometry given in the lower half of the 
figure, As used in the in the Hypersonic Arbitrary-Body Program the skin 
friction surfaces are referred to as an approximate representation of the 
vehicle geometry while, in fact, it has been observed they are as complete 
as generally used throughout the industry for the inviscid calculations. 


From the input element data, the surface normal, area, and area centroid 
coordinates are calculated. In addition, maximum chord length, taper ratio, 
and true area are input for each surface. The latter may be different from 
the calculated area since curvatures have been neglected, The initial surface, 
specified by its maximum chord length and taper ratio, is assumed to be in 
the plane of the skin-friction surface and, therefore, the flow history is only 
approximated. The element planform effect on the average skin friction is 
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b) Representation for Inviscid Calculations 


Figure 10. Geometry Modeling for a Typical High L/D Vehicle 
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included, however, and is discussed separately for larninar and tariestx~ 
flows in‘later sections. The shear force on each surface is assumegiv «. % 
through its centroid in a direciion on the surface parallel to plane contas:..o3 
the surface normal and the free-strea:n velocity vector, as described iu the 
section on free-molecular flow. 


Lecal Flow Conditions 


The required local properties (pressure, temperature, density, and velacity) 
are obtained assuming 2 calorically perfect gas. The pressure on each skin- 
friction surface may be determined from a choice of several of the inviscid 
pressure methods — tangent-wedge, tangent-cone, Newtonian+Prandtl-Meyer, 
and Prandtl-Meyer expansion. At the present time, a continuous strip 
shock-expansion calculation is not available within the Arbitrary-Body Pro- 
gra . and, in this respect, each surfa~e is treated independently of the 
others, 


The skin-friction surfaces and local properties, thus, have been defined 

ina way that reduces the problem of calculating the viscous forces on a 
complex shape to one of solving for the skin friction on a number of constant- 
property flat plates. 


Incompressible Flow 


The basic philosophy behind both the Spalding-Chi and the reference condition 
methods is the same. Namely, that the suitably transformed skin-friction 
coefzicient is given by the constant-property or incompressible formulas 
bas. 4d ona Reynolds number also suitably transformed. To emphasis the 
point, this may be stated ancther way: Th compressible skin-friction is 
given by the incompressible form with appropriate correction factors to 
account for compressibility effects, That is, 


Cr, = Ce/Fy 


) 


Ch = £(Rx;), Rx; = Fry ¢ Rx 


where 
Cg = skin friction coefficient 
Rx = Reynolds number 
( )j = indicates incompressible 
( )s = indicates compressible 
The incompressible formulas used in the Hypersonic Arbitrary-Body 


Program are given in Table I and the compressibility factors, F, and 
Fp, are discussed below. 
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Skin Friction Coefficient, f£(Rxj) 


0, 664// Rx} 1, 328/ V Rx; 


0. 088 (log Rx; - 2.3686) 


Flow Source 


Laminar Blasius 








Turbulent 
(Rxj > Rein} 
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Table I, Incompressible Skin-Friction Coefficient Formulas 


The Sivells and Payne formulas have singularities occurring at low Reynolds 
numbers, However, both occur below the point at which the turbulent values 
cross the respective Blasius laminar curves. Thus, the turbulent incom- 
pressible skin-friction coefficients for Reynolds numbers equal to or less 
than Ryjin are given by the corresponding laminar values. 


Compressible Flow 


Reference Temperature and Reference Enthalpy Method 
Fo = p.jp™ 


1K 1 
FRx = (45/4 Po 
¢c 


where ? is the density, m the viscosity, and the superscript "'*" means 
evaluated at the reference temperature, T*, or reference enthalpy, H”; 








(ee tw Taw 

ae (Al) pe + (A2) T; + (1- Al - A2) 
ok H H 

H WwW AW 

Hy = aes + (A2) H, + (l~ Al- A2) 


The value of the coefficients used are due to Monaghan (Reference 37) for 
Prandtl number equal to 9. 71; 


Al = 0, 5825 
A2 = 0.1875 
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The subscript 'W' indicates the wall vaiue and subscript "AW" refers to 
adiabatic wall conditions given by 





ee eee eee er 
T an = aa 2 ws 
6 6 

where 
Y = ratio of specific heats (= 1, 4) 
M = Mack nurnber 
: \j/n 
r = recovery factor = (P,) 


n = 2 for laminar flow 
n = 3 for turbulent flow 


P,= Prandtl number (= 0. 71) 


Spalding-Chi Method (Reference 38): 





Fo2 A / {arsin(“z) + ARSIN (42) \é 
where te. i Haw | : 

B = = 3M 

C = (A+B)? + “al 


H q Hy, \Pt9 
FRx* (“aw ) |G) |. q= 0.772, p20.702 


Surface Equilibrium Temperature 


In the Arbitrary-Body Program the surface equilibrium temperature 1s de- 
fined as the temperature satisfying the s:2:ady-state heat balance between 
the boundary-layer convection to the surface and the surface radiation to 


sd ea convective heating: QC(T,.) = Ch (Haw - Hwy) 
radiation heating: QR(T,) = Rx te 


where C), is the heat transfer coefficient 
and Ry = €0, €* emissivity {= 0.8) 


o = Stetan-Boitzman constant 
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The surface equilibrium temperature is defined when QC(T,) = QR(TR) 
for Te = Tp. The solution is obtained by a simple linear intercept tech- 
nique illustrated in the sketch and explained briefly as follows. 

Linear relations are assumed for both heating rates 


Qc = AC + (BC)T 


QR = AR + (BR)T 
The four coefficients are initialized as follows. 


1 Assume TCl = TRi = 100 °R 
2. Calculate QCl and QRi 


3. Let QR2 = QC and calculate 


V4 
. (QR2 
TR2 = ( Re 


4, If TR2 > TC2 = Tay, then set 
TR2 = TC2 and calculate new QR2 


QC 1 = QR2 


ee ee eee ree ee eee 


HEATING 
RATE 





TR2 


t C2 = Taw 
TC1=TR} (TC 2 = TR 2, 
(TC 1 = TR, 


TEMPERATURE ———t— 
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The coefficients may now be readily determined and the result of the linear 
solution of the heat balance equation is simply 


T = (AC - AR)/(BR - BC) 


The convective and radiation heating rates are then calculated at this tem- 
perature and checked for convergence: 


1 - QC1/QRi] s EPST, where EPST = 5.0E-4 


If the criteria is not satisfied the cycle is repeated with TC] = TRI = T, 
QR2 = QCl, and TC2 = TR2, The present technique, while lacking sophis- 
tication, is accurate and quite rapid. Normally, two or three cycles are 
required for ideal gas solutions and one additional cycle for real gas cases. 


Real Gas Effects 


It is felt that some comments are in order with regard to the overall pro- 
cedure. Specifically, what is the correctness or justification in using real 
gas reference enthalpy viscous solutions when the local inviscid flow has 
been deterrained only for a calorically perfect or ideal gas? To answer 
this question, an extensive comparison of laminar boundary-layer methods 
was undertaken in support of an earlier study and the details are reported 
in Reference 39. Briefly, the skin friction was determined for the flight 
conditions of the matrix given in Table II, corresponding to the surface 
equilibrium temperatures (emissivity = 0.8) at the one-foot station of a flat 
plate. 


miata Velocity (1000 fps) 
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Table if, Flight Matrix for Skin Friction Calculations 
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Angle-of-attack variation from 0° to 40° in 10° increments and five hound- “5, 
ary-Inyér calculations were made at each condition. The latter correspond ae 
ec -to the combination of three boundary-layer solutions and twe shock wave. 

/ =. _S¢lutions. for losal properties as shown in Table IL 


| Boundary Layer 
Soiution. : 






Exact 





Reference - | 
Enthalpy 
_ - Reference | a teo3 

_. Temperature | — 





2 | 3 





ar ae ree Table HL, Boundary Layer Calculations . 


Alss, additional calculations were made at the flight condition of 20, 000 fps, 
200, 060. feet altitude, and wall temperature equal to 2060°R. 


Methods 1, 2, and 5 are self-consistent with respect to the assump-~ 
tions madé and are reyarded as no.mal calculation modes. Methods 3 

and 4 are inconsistent in the assumptions made between the inviscid and 
viscous solutions and are termed mixed calculation modes, The free-stream 
properties were specified by the 1962 U.S. Standard Atmosphere and 
Sutherland's viscosity formula. The oblique shock-wave solutions are 
accurate to 5-significant digits in the inverse density ratio. For the real 

gas solution, the thermodynamic properties for equilibrium dissociating 

and ionizing air were obtained by the methed in Reference 40. The assumed 
ideal gas is calorically perfect with ratio of specific heats equal to 1, 40. 


‘fhe real gas variation fer the density-viscosity product in the viscous 

solutions was obtained as 2 function of enthalpy and pressure using the 

polynominal equations given in Reference 41, This product is based on 

the most recent thermodynamic data of Hilsennath (Reference 42) and : 

the viscosity calculations of Hansen (Reference 43), The Prandtl nurnber 

was assumed gqual to 0.71 for all the methods. 
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’ : | REAL GAS, EXACT . 
——~——-— | REAL GAS, REF. ENTHALP?® 
mess | IDEAL GAS, REF. TEMPERATURE* 
& ISEAL GAS, REF. ENTHALPY® 
& | REAL GAS, REF. TEMPERATURE" 
*MONAGHAN'S REFERENCE a 
: AT Pp = 07% = 40° a 
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Figure 11. Laminar Skin-Friction Coefficient Comparison 
(Altitude = 200,000 Ft., Velocity = 20,000 fps., Ty = 2000°R) 
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Typical resulis of the comparison are shown in Figure 11. The exact 
scluiion3 were obtained using the Douglas General Laminar Cornpressible 
Boundary-Layer Program as described in Reference 41. The reference 
method calcuiations shown are based cn the coefficient values of Monaghan. 
These were selected 3inte the skin friction calculated consistently gave 
the besi agreement with the exact results. Comparison of the three forrnu- 
lations considered - Monaghan {Reference 37), Michel (Reference 44) and 
Eckert (Reference 45} axe shown in Table IV for the same flight conditions 
#s Figure 11. Major conclusions of the comparison are: 


3, With the exception of possibly zero angle-of~attack the reference 
ternperature method, using existing values for the coefficients Al 
and AZ, is inadequate for predicting ekin friction for the complete 
range of hypersonic flight conditions considered, 


2, The real gas, reference enthalpy methed using Monaghan's forma- 
lation adequately predicts the laminar skin friction over the complete 
flight range considered. The results, however, are consistently 


about 3 to 5 percent lower than the exact calculations, 


3, The mixed calculation mode, ideal gas inviscid ~ real gas reference 
enthalpy is in substantial egreement with the real gas reference en- 
thalpy calculation up to 30° angle-of-attack. 










Reference Angle of Attack in Degrees 


Enthalpy Sond 
0.623|1,956 1.445] 1.753]1.969] 2.056] 2.121 [1.853 1,596 

Micheli 3.24 310,628 1.062;1.447 1.747] 1,53] 2,067] 2,075] 1.798 1.529 
0.243 j0,613}1,038]1.418 |i.77]1 926] 2 042|2,058)1,788 [1.534 


fable LV Comparison of Reference Methods. Values of Cr x 104, 
(Altitude = 200,000 Ft., Velocity = 20.000 fps, Tw= 2000°R) 


























On the basis of the results of this study, the mixed-mode ideal gas inviscid- 
real gas reference enthalpy calculation has been included in the Hypersonic 
Arbitrary-Hody Program. The real gas fluid properties of air are determined 
by the procedures dascribed in detail in Reference 41. Three different formulas 
are used to specify the viscosity. At very low temperatures such as might be 
experienced in a high speed wind tunnel the viscosity is found from the 


68 


Bromley-Wilke results {Refe~ence 46), in the Arbitrary-Body Program these 
are approximated by the fuilowing hinear relationship; 


For T < 225°R 


w= 0, 80383436 Tx107? Wb see 
ft2 


At higher temperatures and for an ideal gas the Sutherland viscosity 
formula is used {Reference 25): 


For T > 225°R 


3/2 
= z x 1078 2bsec 
B= 2,270 Frise. * 10 ne 


For real gas and temperatures greater than about 6000°R Hansen's 
viscosity values are used (Reference 43). 


Viscous ~ Inviscid Interaction 


Under conditions of low Reynolds number and high Mach number, the 
‘rutual interaction of the boundary layer and the inviscid fiow field can 
have a large effect on both the laminar skin friction and surface pressure, 
Boundary-layer displacement effects in hypersonic flow over flat plates 
have been studied at length (e.g., Reference 47) and the present approach 
is limited to consideration of these methods. Basically, a pressure is 
induced from the relatively large outward streamline deflection caused by 
the thick hypersonic boundary layer. The classical approach is to con- 
sider an effective body, made up of the actual body plus the boundary-layer 
displacement thickness, in an iterative solution with the inviscid flow, 

This in itself is an approximation and, in addition, the simplifying assump- 
tions of hypersonic viscous similarity are usually employed. This pro- 
cedure has been adopted for use in the Arbitrary-Body Program and a brief 
background and developrnent of the final equations follow. 


Bertram and Blackstock (Reference 48) presented some simple procedures 
for estimating the boundary layer induced effects on pressure and skin 
friction, These involved the use of hypersonic-similarity~boundary-layer 
theory solutions in an iterative technique with the hypersonic small-disturb- 
ance tangent-wedge pressure equation. The analysis showed good correlation 
with experimental data for surfaces at nearly zero degrees incidence to 

the free-stream. White (Reference 49) extended the theory of Bertrarn and 
Blackstock to include the effect of angle of attack and presented a direct 
method for solving the problem without requiring iterations. White used 
hypersonic small disturbance exp.essions for both compression and ex- 
pansion flows and introduced a new interaction parameter to correlate the 
wall temperature effect. Recently, Bertram (Reference 50) has presented 
more elaborate solutions for the problem employing the techniques of White. 
Implicit to all these solutions is the assumption of a calorically perfect gas 
and a Prandtl number of unity. 
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White's solution has been uzed in the present analysis because of the rela- 
tive simplicity in its application, His numerical results showed the local 
pressure to be nearly a linear function of the iteraction parameter, }; 


P= Py (i+ B) 


where d 
B= m Pp ’ 
© 
and 3 u/2 
: GM ( c } 


Jivay \R® 


The quantity G is a simple function of wall temperature and specific heat, 
C is the Chapman~Rubesin viscosity coefficient, and j is the Manger 
transformation parameter: two-dimensional flow, j= 9; axially-symmetric 
flow, j=l . 


W 


In the above equations, P is the local pressure to tree-stream pressure 
ratio, and the subscript "o" refers to the inviscid value obtained frozn the 
hypersonic small-disturbance relations. 


Bertram's (Reference 48) correlation far local skin friction coefficient is 


1/2 
Cy = 0,664 Ky (FS) 


where Kj is a pressure gradient and wall temperature correction factor, 
The shear on the surface is 


Ty = fase dA 


In the present analysis, the approach taken is to determine the effect or 
factor due to viscous-interaction using White's method and then to modify 
the previous result without interaction accordingly, This viscous-inter- 
action factor, Kyy, is obtained by carrying out the integration of the pre- 
ceding equation and is defined as follows; 


(Twiyt ——} Vl+ Bey +1 
Ty = Ji+ Bey + Bo, log, Ba age 
Cr 
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where Be, is based on the root-chord and K, has been assumed equal to 
one, This expression is for a plate with taper ratio one, but the integration 
could have been done for an arbitrary value (e.g., Reference 51). In the 
present application the planform effects are included in the shear force 
without interaction, Ty. This application results in a slightly lower factor 
but has the advantage of permitting a step-by-step build up and comparison 
of the overall viscous forces. The magnitude of the skin-friction correction 
facter using the above techniques is shown in Figure 12, 


The induced pressure on a surface is determined as an increment in pres- 


sure coefficient, 


Pp 
AC 2 C=C) iS 
Pp P Po 3 & 


The average pressure increment, P - P,, is found by summing the local 
pressure distribution over the surface. 


vs vd 
B-p- + fip-p)aa 
Substituting the expression for local pressure and integrating gives 


Bp - B = 2émacy 


The AC, due to induced pressure is determined for the skin-friction geo- 
metry répresentation of the vehicle shape and effects due to the planform 
shape and due to the initial surface are discussed in the next section, 


The basic hypersonic smali-disturbance relations for calculating pressure 
are: 
For compression flow (K 2 0) 
/ 2 1/2 
2 Y+l 2 , ¥+1 
Ped fT) i + nc i+ (HE x) 
For expansion flow (- 2/(y- 1)s Ks 0) 
27/(¥ -i) 
P= [: + 1 K] 


The similarity parameter, K, is given by; 
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Figure 12. 
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HYPERSONIC INTERACTION PARAMEYER ~ X 


Effect of Viscous Interaction on Skin Friction Coefficient 
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where K, = M,siné (§ is the surface impact angle} and Ky, a boundary- 
layer growth parameter, is taken equal to 1.0. 


White (Reference 49) observed that the pressure equation (either compres- 
sion or expansion) and the expression for K constituted a first-order 
nonlinear differential equation in P(,) and obtained numerical solutions 
directly without iteration. The results are shown in Table V from which 
White also observed that the pressure could be approximated by the Jinear 
relationship 


P= Pp, + md 


where P, and the slope parameter, m, are just functions of Ky. FP is 
given by the hypersonic similarity relations as a function of Ky and, in 
the Arbitrary-Body Program, m is approximated to the data of Table V 
by the following analytical curves: 


For -2/(y-1)s K < -3.9, 


m = 1,424 + 0.219 K, 


For K 2 -3.0, 


m = 1.9156 + 0.41727K, - 0,0419101K," 


- 0.010427K,° + 0,00214381Ky* ~ 0,000103217K, 


5 








Similarity Parameter K, 


473| 8.734 | 44.14 | 170. 


~ 
i] 
Ww 








Q.0 | 0.002 | 0.028 | 0.210] 1.000] 3. 2 
0.5 | 0.173 | 0.339 | 0.748} 1.835] 4.555] 9.930} 45.41 | 171.4 
1.0 | 0.428 | 0.736 | 1.379} 2.777 | 5.722} 11.18 46.70 | 172.7 
1.5 | 0.738 | 1.192 | 2.059| 3.740 | 6.914] 12.47 48.01 | 174.0 
2.0 | 1.092 | 1.695 | 2.770) 4.709 | 8.108] 13.76 49.33 | 175.3 
2.5 | 1.485 | 2.234 | 3.506} 5.679 | 9.294] 15.07 50.66 | 176.6 
3.0 | 1.908 | 2,801 | 4.260; 6.651 }10.47 | 16.37 Seek mae a a re] 
32 9),| C2399 | 3.592: | 5029-17. 622-401. 64 1 47.67 53.34 | 179.3 
4.0 | 2.833 | 4.004 | 5.810; 8.593 112.80 | 18.96 54.70 | 180.6 
4.5 | 3.328 | 4.632 | 6.601! 9.505 413.95 | 20.25 56.06 | 181.9 
| 5.0 3. 840 5.275 | 7.400] 10.54 115.09 | 21.52 51342. | 183.2 


Table V. Numerical Solutions for Pressure Ratio P (y =1-4) 
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Planform Effects 


The previous sections have dealt with the determination of the local skin- 
..friction coefficient or the average skin-friction coefficient per unit span, 

In this section, the determination of the viscous force contribution of a 
surface element having a planforrn shape of the type shown in the sketch 
below is considered. In the derivations that follow it is implicitly assumed 
that the root and tip chords are parallel te the oncoming flow. 


FLOW 
DIRECTION 





The product of local skin-triction coefficient (Cs,) and dynamic pressure 
(4) is integrated over the surface area (A) to obtain the shear force: 


Ww fas ME ae 


{The symbol T is customarily used to define shear stress, however in 
the present text it is used consistently as a force. This is done to avoid 
the unnecessary use of area ratios in the defining equations and at the same 
time retain the significant connotation associated with the symbol.) 


T 


The shear farce cn: ach surface is then written as a coefficient with re- 
spect to the free-stream dynamic pressure (q.) and a speci:ied reference 
area (S), T. 
W 
Cyr = 


fo] 


n)- 


and summed over all surfaces to obtain the vehicle characteristics due to 
viscous forces. 
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Laminar Shear Force 


The local properties ace constant on each surface and the abave ecprea sian 
becornes 


- 
~ 


) 
b = 
é 
Ty = 4 (Cede eee (x) re. (a 
a 


S 


where the surface has root chord c,, span b, and (C; Jey is evaluated at 


the root chord. The local chord length may be expre ese as 


c 2c, (1-- TR)o] 


where TR is the taper ratio (= ct/c,;) and 9 is the normalized span 
dimension (= y/b). Substituting this expression and cornpleting the inte- 
gration gives the shear force on the surface as 


Al 1+ TR+ VTR 


Tw = gq, A(Cr,) = 
wey ese Fo'cr 313 TRU + JTR) 


where (Cr, Pp is the local, length-averaged skin-friction coefficient 


evaluated ake root-chord. 


In the Arbitrary-Body Program the shear force is expressed in terms of 
an average chord length, Cc; 


where 2 


ae [teres TR 
ae (1+ TR)(1 + /TR) 


(eo) 
t 


Viscous-Interaction 


As was explained in the previous se&tion, the effect of planform on the 
shear force is not determined directly for flows ...th viscous-interaction 
but is included in the calculation of shear force without interaction. This 
procedure results in a slightly lower force but has the advantage of per- 
mitting a step-by-step build-up and comparison of the overall viscous 
forces. There is, however, an additional effect on the induced pressure 
due to planform shape which 1s accounted for, 


The average pressure is obtained by integrating the local pressure over 
the su-face: 
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where — = x/c,, the normalized streamwise coordinate, 


Substituting the expressions for 


cyb 
A= eo (l+ TR) 


c/e, =1- (l- TR)y 


the integration is easily completed. The result is 
Py = Po 1+ 9 Be. 1+ TR+ VTR 
(+ TR)(L+ (TR) 


where 


The average pressure increment for the surface is then 


P, - Pot : mc | _1+ TR+ VIR 
(1 + TR)(1+ TR) 


which for TR =1 reduces to the value previously given. 


Turbulent Shear Force 


Because of the nature of the assumed skin-friction formulas, a different 
approach than used for laminar flow is taken to obtain the turbulent shear 
force. The end result, however, is an approximate solution which is very 
similar to the laminar result. The shea+ force equation is derived as 
follows. 
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Ty = fa, Cy, dA 
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The variable of integration is transformed to the local chord-lengfh Reynolds 
number in two steps, First in terms of the chord length c, 


Cr c Cr, 
Tw = qs b obs TRY: = de 


Cy 


Next, the variable of integration is transformed to the incompressible 


r : ‘ 
Reynolds number, Rc; = Fry (782), and normalized with respect to root- 
6 


chord values; 
Cr : 
woe MBS Fs) 6) oe 
Wo - TR) de Vow MOF eriy \ROrhs 





Cyb 
Noting that the surface area is A= “— (1+ TR), and also tat aoe) = 
F ‘ Cy 3 
ce). the shear equation becomes 
ZL 


“r 
1 


2 Rc Cr Re 
T = A(Cc ees f peta d 
W qs ( Fs Ne (; ee me) (FE) (<=) (a) 
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With a simple power-law skin-friction formula this equation is easily 
evaluated; 





1 
Cc “WN 
(c= = (Rs) , Where N is positive 
and 1 1- 1 
Ty, = A(Cr.) c d/ Rc 
wt 4 Fy)e,( z) f (Re.\ ( 
r\l-TR TR \Rey/ Re, }, 
1. TR2- 
= q A(Cpr ) - \||—————— aes 
5 bor I 1-TR 


For laminar flow N= 2 and it is easily verified that this expression 15 
identical to the one previously presented, 


in general, the skin-friction coefficient is not given by a simple power-law 
relationship and this is the reason for deriving the turbulent shear with the 
Reynolds number as the independent variable. 


The use of the Sivells and Payne formula in the shear equation introduces 

a singlularity in the integrand and the function is nonintegrable. However, 
this singularity occurs at a Reynolds number much below the laminar cutoff 
and the shear equation may be integrated numerically. Several examples 
for the numerically determined integrand are shown in Figure 13. The 
upper-bound ‘represented by laminar flow and a lower-bound represented by 
constant skin-friction are also shown. The curves are smooth and the area 
under each curve times the quantity Z/(1 - TR“) is the factor by which'the 
shear increases due to a tapered planform. 


It may be observed from Figure 13, that even with a large variation of 
Reynolds number on the planform (for example, Rc, = 10 9 to zero at the 
tip), the major contribution to the integral is obtained over the first decade 
(Rce/Rcey = 1.0 to 0.1), In the case of the upper-bound (laminar flow) and the 
lower-bound (N =@) this contribution is 97 and 99 percent, respectively. 
This then, suggested the approximate approach of representing the Sivells 
and Payne formula in the integrand over the entire Reynolds number range 
by a local power-law fit obtained as the average over the first decade. 


Thus, the shear on the surface is obtained from the power-law solution 
with the exponent parameter, N, given as (for Sivelis and Payne); 


_ log Rey - 2 
~ — 0,8686 


Alternately, as was done for laminar flow, the shear force may be expressed 
in terms of an average chord, c; 


*w = 45 A(CEy). 


where 


and 
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Figure 13. Planform Effect on Shear Force 
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Initial Surface Correction to Shear Force 


When an initial surface is specified, the shear force is determined for the 
combined surface geometry, for the initial surface, and the difference 
obtained as the value for the surface of interest. This in effect is dealing 
with three surfaces which have the following characteristics (see sketch 
below): 


1, Initial surface; Area Ay, maximum chord length Lj, taper ratio 
TR), and shear force Ty). 


2. Surface of interest; Area Az, maximum chord length L2, taper 
ratio TR2, and shear force TWo> 


3. Combined surface; Area A3 = Aj + Ag, maximum chord length L3, 
taper ratio TR3, and shear force Tyw3. 


The shear force on surface 2 is 


TWo = Tw 7 TW 


= (Cr Kyy) 
= dg Aa(Cr, Kyy); 41 - rks ~ 


In the Arbitrary-Body Program this is compacted to the form 
TW = qs Ag (CR, Ky1)3 (l- FF) 


where FF has the mnemonic form factor or friction factor. Three possi- 
bilities are considered in determining the friction factor: (1) both surfaces 
laminar, (2) first surface laminar and second surface turbulent, and (3) 
both surfaces turbulent. 






/ 
7 _ INITIAL SURFACE 
ra AREA A, 


SURFACE OF INTEREST 
AREA Ay 


ania cee ian 
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Initial Surface Correction to Induced Pressure 


The average pressure on surface 2 is defined as follows: 


where F; is the force on surface i. The average pressures on the initial 
surface and on the combined su: face are given by 


1+ TR) + ¥TR | 


m2 °P hig coe, pennies cas cat 
’ a 3 Fe TR) + YTR}) 


~ 


1+ TR.+ /TR 
ee ee re 
(L+ TR3)(1+/TR3) 


and the areas by 


Ay, = bly (1+ TR,)/2 
A> = bLa(i+ TR2)/2 
A3 = bLa(l + TR3)/2 


Substituting these expressions into the above definition and after same 
algebraic manipulation the result may be writien as 


l 
B.-p -$oa(e)} 1+TR, +¥TR, \f (2) (ee =)| 
zd re) - PaaaEDennCanay (puunnieeeninds 

2/1 UeTRa+VTRs)] | \E3/ TR TRALVTR, 





The length L, is defined as the maximum chord length of the combined sur- 
face, soas L)~-®0 it is readily verified that the pressure reduces to the 
same expression previcusly given for a single, tapered plate. 


Viscous Force on Blunt Bodies 


The earliest space capsules were designed with large spherical nose 
caps and flew ballistically at zero degrees angle of attack. For such 
vehicles, it was found that inviscid flow field calculations were ade- 
quate to p* sdict the splash point. The later generation capsules were 
designed to fly at angle of attack to provide lift and it has been shown 
that viscous forces can have a significant effect on predicting the 
splash point. The theoretical solution, then, must provide some means 
for estimating the viscous effect. 
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The procedure included in the Arbitrary~Body Program is that developed 
by Goldberg of the General Electric Company (References 52 and 53), 
This method is given in the form of relatively simple correlation formulas 
in terms of the shock-layer Reynolds number and inverse density ratio, 
The method is applicable to the low density conditions associated with high 
altitude entry and is equally suited to real gas or ideal gas analysis, 


The shear force in the stagnation region of a blunt-faced body is given as 
Tw = TWo Kyy 


where the shear without low density or viscous-interation effects is 


Two = ; IoA 2 cos 6 
(1 - 0,455 fe) Le 
and 6 is the surface impact angle, 
€ is the inverse density ratio, = (?2/A, ‘ 


Re is the shock Reynolds number 


P2 U2Rp/ U2 


Rp __ is the body nose radius. 


The vigcous~-interaction correction factor, Kyy = Tw/ TWo: was obtained 
from higher-order analysis of the boundary-layer flow (Reference 52 ). 
The present authors have developed a correlation formula to represent 
these solutions in the Arbitrary-Body Program. This factor, a compli- 
cated function of both shock Reynolds number and density ratio, has been 
approximated by a combination of exponential transition functions of the 
type described by Grabau (Reference 54). These are 


ae 1 
even transition: y = “Sarcxp Rie ae 
~ to 
ake 1 
odd transition: y = Trop kOe 
° 


These functions are essentially the kernels for the Bose-Einstein and for 
the Fermi-Dirac distribution functions, respectively, for the even and the 
odd transitions. The notation of transition is used since these functions 
represent the smooth transition from one asymptote to another; the even 
case does not have a point of inflection and the odd transition has a point 
of inflection. 
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In the present application, a correlation formula for the viscous -~inter- 
action parameter has been obtained by a combination of an even and odd 
transition function. The curve is considered to have three asymptotes 
(see the sketch below); Yi, Fl, and F2. First an even transition is 
determined for the curve between Fl and F2 and this is designated Y2. 
Next, an odd transition is established between ¥l and Y2. The curves 
are adjusted through.the values specified for the exponential constants, 
K, and the origin coordinates, X 9. Details of this procedure are given 
in Reference 54. 





The correlation formulas developed for the present case are as follows. 


Independent variable X = log (€3Re,) 


Fl = Al + Bl(X) 
Al = 0.667 
Bl = 11111 

F2 = 1.0 

Yl = 0.0 
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Y2 = Fl + COTS EVE Ie KGET) att 
T0- exp [EVK (K - XOEV) 


EVK = -1.80 


1.0 + exp [ODK (X - XOOD)] 


ODK = -2Z.0 
XOOD = AOD + BOD(log e€) 
AOD = 1.0 


BOD 


3.2907 


Comparison of this correlation and the boundary-layer solutions are 
shown in Figure 14, The general shape of the curves is well repre- 
sented by the correlation, although some accuracy is lost, particularily 
at the peak of the €= 0.04 curve. It would be possible to tailor-fit 
each of the €-curves through further variation in Fl, the exponential 
constants and origin coordinates. However, since only three solutions 
were available, the determination of more accurate fits was not deemed 
justified. Three additional €-curves are given on the figure to demon- 
strate the behavior of the correlation formula. , 


An example of this technique is shown in Figure 15 were the predicted 
values of lift coefficient for the Gemini space capsule are compared 

with experimental results (Reference 55), The modified Newtonian 
calculation has been performed for the entire shape and the viscous 
calculations (broken lines) made only for the blunt face. The present 
comparison, due to the limited data used, may not completely justify 

the method, but it does show the significance of the viscous contributions. 


The blunt-body viscous calculations are not limited to entry capsules 
but may be applied to any blunt portions of a vehicle (e.g., leading 
edges), The method is primarily dependent on impact angle and, there- 
fore, the detailed inviscid geometry is used. It is for this reason that 
the method has been included as one of the inviscid force options. Zero 
contribution is assumed for shadow flow. 
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Figure 15. Gemini Lift Coefficient Comparison 
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Control Surface Forces 


An important feature of the flow about hypersonic control surfaces is the 
boundary-layer flow separation phenomenon. Flow separation on the control 
surface and on the surface of the vehicle ahead of the control can have a pro- 
nounced influence on control effectiveness. This is a very difficult problem 
to analyze theoretically with any degree of accuracy. However, the use ofa 
simplified flew model and empirical boundary-layer separation data will 
allow the solution of this problem with sufficient accuracy for most prelimi- 
nary design purposes. 


The basic factors involved in the control surface flow separation phenemena 
are illustrated in Figure 16. In these flow models the pressure in the 
separated region is taken to be the plateau pressure. This region is well 
defined for laminar flow and may affect a large portion of the vehicle surface 
on, and ahead of, the control flap. For turbulent flow the separated region 
is much smaller but, although a true pressure plateau may not exist, a flow 
model using the inflection pressure will produce useful results for most 
purposes. 


The boundary layer separation at point Xgpp is termed a free-separation 


and the separation angle in this flow model depends only on the local 
Mach number (Mj) and Reynolds number (Ry) ahead of the separation re- 


gion. The extent of the separation region as indicated by the separation 
length Lorp is determined by the strength of the shock at the flow reattach- 


ment point Xp. As the control surface deflection angle ¢ is increased, 
the flow turning angle at the reattachment point increases, the final overall 
pressure rise after the reattachment shock increases, more flow is forced 
back up the boundary layer and into the separation region, and the length of 
the separation region increases to accommodate the greater flow in the 
separated area. The length of the separation region continues to increase 
until the reattachment point reaches the control surface trailing edge. 


Present theoretical methods are completely inadequate for the prediction of 
separation effects; complete reliance, therefore. must be placed on empirical 
correlations. The relationships used in this program were taken from the 
work of Popinski and Ehrlich in Reference 56. The basic approach used in 
the Arbitrary-Body Program is essentially the same as this reference. 
There are, however, a few significant differences in the details of the pro- 
gram application. Since the program can perform shock-expansion calcu- 
lations this capability has been exploited in the flow-separation application. 
Also, since all of the calculations are carried out by the computer program 
rather than being done by hand it is possible to use a strip theory approach 
that permits a rough assessment of the three-dimensional effects, The 
computations within the program assume an ideal gas. 


The application of the flow-separation criteria in this program also has 
another very signifi-ant different from Reference 56. The program uses 
shock-expansion theory to determine the flow separation pressure changes 
and then applies th:se separation corrections to the pressures calculated 
for the vehicle, usmg any of the other force methods. This process is 
illustrated in Figure 17. 
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Figure 17, Correction of Normal Pressure Method 
Results for Separation Effects 


The final surface pressure is calculated by the relationship shown below. 
AP 


6 
Pp = P. + s P 
N Pen ( , 


This technique not only gives a smooth trend of data as the deflection angle 
goes to zero, but also insures that 1t will match the basic vehicle results 
using whatever normal force method has been selected. 





The essential features of the flow-separation calculations are given below. 
The notion is maintained consistent with Reference 56, 


The calculation of flow-separation effects uses a three-cycle process, On 
the first cycle a streamwise strip of elements is analyzed using the shock~ 
expansion method to determine local conditions at the hinge point. The 
inviscid pressure rise onto the flap is calculated also and compared with 
the incipient separation criteria. The transition Reynulds number is input 
to the program. 


On the next cycle the shock-exparsion method is again applied to the same 
streamwise strip of elements, only this time sevaration effects are ac- 
counted for. This is accomplished by checking at each element to see 
if the flow separation point has been reached, and by calculating the 
new pressures. Again, the equations in Reference 56 were used in 
these computations. For other reviews of the separation problem see 
References 57, 58, and 59, 


On the third cycle of calculations along the streamwise strip of elements 
the pressures are calculated using the normal input pressure calculation 
method (IMPACT and ISHAD). At the same time, these pressures are 

corrected by the flap deflection and flow separation pressure increments. 
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The above discussion gave a brief outline of the techniques involved in the 
flow separation calculations, Of course, a thorough study of the listing of 
the FLOSEP subroutine is necessary to obtain a thorough understanding of 
the methods used. The following description of the equations used will 
assist in this study. , 

The basic purpose of the first flow separation analysis cycle is to dete -mine 
the local flow properties at the hinge-line element, and with this information, 
to determine the pressure rise in going from the fore-surface element at the 
hinge line to the first element on the flap. This pressure rise is then com- 
pared with the empirical incipient pressure rise equation to determine if the 
flow separates, 


On the first FLOSEP cycle the flow properties on each element in the stream- 
wise strip are calculated by the shock-expansion method, This gives complete 
freedom to analyze curved vehicle surfaces in front of a flap (the fure-surface) 
and also on the flap itself. In many such cases this approach will give useful 
results even though the empirical separation criteria equations used are 
derived from tests on flat surfaces only. 


When the hinge-line element is reached on the first FLOSEP cycle the program 
checks for flow separation. The flow turning angle at the hinge line caused by 
the control deflection is calculated from the following relationship. 







2 2 2 
=, /(N.N -N,Ny )°+ (NUN, -NSN, )° + (N,N -N,N. ) 
{ Y Zr Zz Yay zZ Sy Xx Zu x Yu Y Xe 
where ‘ 

N., N.,, N = the direction cosines of the first element on 

x’ Y Zz 

the flap. 

Ny, , N, , N = the direction cosines of the hinge-line element. 

xy Yu “ey 


Note that when the flap is not a plane surface the flow deflection angle cal- 
culated by the above equation will be different than the program input control 
surface deflection angle. 


The inviscid pressure rise at the flap is calculated by the shock~expansion 
method (oblique shock or Prandtl-Meyer expansion as the case may be), 
This pressure rise is then compared with the incipient pressure rise cri- 
teria given in Reference 56. 


2 -0. 306 
2.93 (Mq_ -1* 
(Cy } = (R 1/4 Laminar 
a@ inc ey 
HL 
2.2 
(CL) = 1/10 Turbulent 
Py inc Re ) 
OL 
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where 


(Cy ) = the pressure rise required to cause incipient 
@ inc separation. 
My = the local Mach number at the hinge line. 
R = the Reynolds number at the hinge line based on 
te, the local hinge-line flow properties and the 


distance from the leading edge. 


Subsequent FLOSEP calculation cycles require a value for the flap chord 
length (Crap) This length is taken as being the distance from the hinge 


line to the flap trailing edge as given by the following equation. 


2 


3 
yo + (Zp > 27) 


Gi, Ke ae 


flap seas 6 is 


(Yor - Yury 


where 
Y 


Zone = the average values of the X, Y, Z element 


XqE 
coordinates at the trailing edge. 


TE’ 
Viscous separation effects are calculated on the second FLOSEP calculation 
cycle, The general procedure used differs from the method of Reference 56 
in that no iteration process is used to determine the separation point. Since 
the vehicle surface is not restricted to flat surfaces this technique is not 
possible. Instead, as the local flow properties on each stream-wise element 
are calculated a check is made with the appropriate equations to see if the 
separation point has been reached. When flow separation 1s indicated ona 
given element a linear interpolation of parameters between the last non- 
separated element and the separated-element is used to determine the exact 
point and conditions at the separation. 


The following calculations are required to provide information for the check 
for flow separation. The local Reynolds number at the element centroid 
is given by 


R = X (R_ / ft) 
Cox! LE Ly 
fr) 
where 
XUR = the distance from the leading edge to the element 
7 centroid, 
R_ / ft = the Reynolds number per foot based on the local 
e ony 
a element flow conditions. 
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The plateau pressure rise due to separation is given by the following 
equations. 


1. 86 (M,/ . 17? 262 


(Cy ) = . Laminar flow 
a Pp 1/4 
RR) 
ox ! 
,o) 
-0. 309 
(Cy ) 1.91 (M, 1} Turbulent flow 
a Pp 0,1 
Re) 
ax! 
° 


The ratio of the plateau pressure to the local element pressure upstream 
of the separation is given by 


P_/P = 0.7M“ (C 
oS a 


p + 1,0 


The ratio of the distance from the start of separation to the hinge line (dy ), 
and the local boundary layer thickness (§ o> is given by 


iu 


5 -4, 1 


5.69 X10" M, 3.5 


d,/ § (Pp/P, - 1) Laminar flow 


ce) 


1.1 x 10° Mee (P/P, - 1) Orne” Shaaban’ How 


Uy 


dj og 


The calculation of the boundary layer thickness requires the local surface 
Reynolds number per unit length based on the reference condition (Rey / ft). 
This parameter is obtained by the temperature subroutine (TEMP) used in 
the skin friction calculations. The surface wall temperature required to 
find the Reynolds number may either be ietermined by the temperature 
routine or input. The same options used in the viscous force calculations 
are available to FLOSEP. 


The boundary layer thickness is given by the following equations, 





Laminar 
86 we se 
a = (R * rey ll Turbulent 
a 
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From the above equations for d,/é, and 6o/Xur we may obtain the distance 
from the start of separation to the hinge line (d,). The separation distance 


parameters are shown in Figure 18 as taken directly from Refererce 56, 


PRESSURE DISTRIBUTION a7 REFERENCE LINE 
ALONG THE SURFACE : Pt Bees oe 
Do Pp fo 
Ps 
P2 


DISTANCE 
ALONG THE 
SURFACE 





Figure 18. Definition of Interaction Parameters for Separated Flow 


The check to see if the flow separation point has been reached is accomplished 
as follows, 


Xspp = X -d 


AXcorp = X -~ X 


If Xorg = Xopp: then the separation point has been reached. 
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where 


x 


LE the distance from the leading edge to the hinge line. 
H 


x the distance from the leading edge to the separation 


point. 


SEP 


The exact flow separation point and the attendant flow properties are cal- 
culated by a linear interpolation between the data on the element before the 
separation and the actual separation element. 


The procedure used in calculating the other required flow separation geomet- 
ric parameters and the separation pressures is outlined below. 


The ratio of the free interaction length (255) to the boundary layer thickness 


( &,) is obtained from the following equation. 





fee 
“fis 2.47x 10° M.“*? (Pi/p -1)?*? ® Lamina: 
re a P'"o 
Le: ~8,4 
Sond, sca 1,84 x 104 Pp/P, =o Turbulent 
bo mee 

(og 


The free interaction length is now calculated, 


hes 


i = (Ng /66) 8 


If the flow is separated on the first element of a strip the parameter Bes is 


set equal to 0.0. Also, the parameter X is set equal to 0.0. 
LE sep 


The ratio of the downstream interaction length to peak pressure (¢,), to the 


upstream interaction length (d,) is given by the following equations. 


For Laminar Flow 


Cra ds 
= >i1: an = ,545 - 04 (M 6») 
1 1 
Cra a5 
tap < 0.25 : -— = .273 - .02(M_5$,) 
d d af 
1 ] 
M6. = 5: S2/4) 2 , 344 
af 
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where 


Cha Cas 
=, if .25 < or hae sl 
i 1 


_flap_ ; = : 
> J: d,/d, 1,16 0. 33 M5 ¢ 


fla - ; 
FP < 0,25: dj/d, = 0,58- 0.165 M8, 


Bae Ue 37 





where 
Cra 
1, for —~P. z 1 
1 
Ona fla 
Ts —7=., for 1 2 ~y~2 2 0,25 
1 ] 


Cra 
0.25, for 2 <s 0.25 
1 


From the above equations the downstream interaction length to the peak 
pressure (d,) may now be calculated. 


The angles associated with the separation are given by 





@ = sin7! (Cy ) oe + —, 
ap M 4 


2 2.0 sin 6 
tan > = / (Gor - 9) moe e | 
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where 


e 


the flow deflection caused by separation, 
@ = shock angle for the flow turning angle, $. 
If the flow has separated at the leading element the angles are given by 
d, sin § 
-1 | "2 f 
> = tan oO 
qd, + d, cos 5, 


The ratio of the downstream interaction length to pressure rise (d to the 


), 
upstream interaction length (44) is given by 2 


d, 7 tan > 1 

dy ~\tan bp - tang COS 
and 

d, = (d,/d,) dy 


If the distance qd, is greater than the flap chord length, Chap’ the separated 


region as determined fron. the above equation extends beyond the end of the 
flap. When this condition occurs the distance d, is set equal to the flap chord, 


Crap’ and a new value for qd, calculated based on conditions that existed at 


the originally calculated separation point. 


ie cae as 
This value for qd, will be smaller than the original value and wiil give a new 


flow geometry model with the same flow turning angle, 4, but with the 
separated region reduced to extend to the exact flap trailing edge point. 


The plateau pressure caused by the separation turning angle 9, (Cy ), is 
calculated by the oblique shock compression subroutine. Ps 


The flow turning angle in going from the plateau pressure region to the final 
peak flap pressure area is found from 


$, Shy ae 


The resulting peak flap pressure (Cy ) is also calculated by the oblique 
shock compression subroutine. = 
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The flow separation pressures are distributed on the fore-surface and the 
flap by using the following equations. 


Lt < X C = C 
Be TP eee PsEP * 1, 
If X > X and < (X + £..) 
LE LEgep LEcpp fi 
“Py. Px “Py. Px 
Cpe 7 (| —2-——_] “ie * Sp, “| —4——- “Le 
‘ Lf Bei 
where 
Cp = the surface pressure at the separation point. 
Xx 


The ivllowing equations are used to distribute the separation pressures on 
the flap. 


If X < (X + d.) Cc =C 
LE LE 2 Popp Py 
P 
If Xur > (Xin + d,) and < (Xp + d,) 
H H 
°P, : “P, 
a (ae BP) Gia + Sige Sida CS 
SEP d d l 
a 83 P 
For d,>C Cc = °C 
2 flap Pots 20 
If d, >C, Cc ae ©, 
3 flap Popp es 
If X,, 2 (X + d.) C = C 
LE LE, 2 aa a 


At the end of the second tlow separation cycle the above pressure coefficients 


are converted to pressure with the following equation. 


2 
P = lc VoM< ot 10)? 
eeuie: || SER 2 © 
where 
Po => the free stream pressure. 
M, = the free stream Mach number, 
Pe = the pressure on the surface using shock-expansion 
SE : p : ‘ 
vis methods and including viscous separation effects. 
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The purpose of the third cycie of calculations down a strip of elements is to 
determine the basic surface pressures using the normal program pressure 
calculation method (i.e., Newtonian, tangent-cone, etc.), and to determine 
the final surface pressures including the viscous separation effects, For 
this cycle the control surface is set to the undeflected position. The final 
surface pressure is given by 





4P.. 
P= Py 4 7 e Pyy 
SE 
where 
2 = the final surface pressure including contro! surface 
effects, 
Pr = the surface pressure calculated by the normal program 
input pressure calculation method. 
AP, = the change in surface pressure due to control surface 
e deflection and including separation effects, 
Por = the surface pressure using shock-expansion methods 


without the separation effects. 


Typical results obtained with the above analysis techniques are shown in 
Figure 19, 
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Effect of Flow Separation on Surface Pressure 


Figure 19, 
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Propulsion Effects 


The operation and design of an air-breathing propuision system ofa 
hypersonic vehicle can have a strong zafluence on the vehicle stability 
and control characteristics. For some aircraft, such as scramjet- 
powered vehicles, the engine, including the inlet ard exhaust systems, 
may be highly integrated into the vehicle design, For these configura- 
tions the bcokkeeping system as to what is engine and what is airplane 
is no longer easy to resolve, A cetailed analysis of these problem 
areas is obviously beyond the scope of this program. These problems 
must be solved by detailed ergine-airframe studies. 


It is important, however, that the Arbitrary- Body Program be capable 
ef properly using the results from these propulsion system studies in 
evaluating a total vehicle's stability characteristics. The approach used 


‘is outlined in the following discussion. 


In investigating the propulsion effects on hypersonic stability it is 
usually desirable to examine each of the various components of the 
engine thrust system. The equations for a general air-breathing pro- 
pulsion system are given below. 


T e Pe (V2) - m V2 bi (P, 7 Py) Ay 


R = ees (Vv) 5 (Pp, = P,) Ay 


where 
= grogs thrust 
R = ram drag 
m, = mass flow of air through the engine system 
TH, = mass flow of fuel added to the stream 
Vy = flow velocity at the forward flow control boundary 
Vv, = flow velocity at the rearward (exhaust) flow control boundary 
Po = reference pressure 
Py = pressure at control boundary i 
P> = pressure at control boundary 2 
A = area 
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The vector relationships for a simple engine system are illustrated 
in the diagram below 


Gross Thrust 





The propulsion system has two rnain influences on a vehicle's stability 
characteristics. First, the engine ram drag and gross thrust vectors 
create moments about the vehicle center of gravity. Second, the ex- 
haust jet may at times expand onto the aft surfaces of the vehicle and 
alter the local pressure distributions. The numerical values of the 
parameters involved depend upon the engine-airframe bookkeeping sys- 
tem adopted. The final result, however, should be the same under any 
system, Any given bookkeeping system serves only to insure that all 
factors are properly included and that no factors are included twice. 


The engine related forces are input to the program in the form of vectors 
positioned in space. The input information includes the magnitude of the 
force in pounds, its point of action on the vehicle (in the vehicle X, Y¥, Z 
coordinate system), and its direction. The direction is given by the force 
vector direction cosines. Any number of vectors may be used. 





The force vectors are converted into vehicle cvefficients by the following 
equations. 


AC, = F(-NS)/(a5,,. 0) 


ACY = F(-Ny)/(qS_ 0) 

AC, FN, )/(aS. 0) 

ACg ACV(Z 0 ort - Zggl/ SPAN 
POC HY ont 7 @o,)/SPAN 
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AC, = AG AX sont - Mag (MAC 
FAC A(Z ont - Bagi MAC 
AC, = AC U(X ont « Regi ee N 
~AC AY ont - Yel PAN 
where 


F = force vector in pounds 


Ny, Ny; Ny are the direction cosines of the force vector 


q = free stream dynamic pressure 


These coefficients are converted to the lift and drag directions by the 
standard transformation equations previously discussed. 


Dynamic Stability Derivatives 

The dynamic stability derivatives due to pitching velocity (C ) and 

vertical acceleration (c,.) are important in the damping of the short- 
Qa 


period stability mode. For flight at hypersonic speeds the derivative 
Cus provides the major part of the damping since Cas _ is usually 
a 


quite small. 
The derivative Cu is the change in pitching moment coefficient with 


varying pitch velocity and is commonly referred to as the pitch damping 
derivative. This derivative represents the effect of rotation of the 
vehicle about a spanwise axis at constant angle of attack such as ina 
steady pull-up. Because the vehicle is rotating, different parts of the 
vehicle have different velocities relative to the free-stream depending 
upon their distance from the center of rotation. This is essentially a 
steady-state problem and may be solved by the use of steady flow con- 
cepts. The common definition for Con is 

q 


(=) 


a 
m 


™a ay) 
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where 


Co. = pitching moment coefficient 
q = pitching rate 
é = reference cherd for moment coefficients 


(frequently the mean aerodynamic chord) 
Vv = flight velocity 


The method used in the Arbitrary Body Program to calculate this 
derivative has been presented earlier in this report and will not be 
repeated again in this section. The basic approach used involves the 
calculation of the local relative flow velocity for each part of the ve- 
hicle depending upon its distance from the center of rotation and the 
rotation rate. Once the local velocity direction is known the surface 
pressure is calculated using any of the methods available in the progra. . 
as desired (except the shock-expansion method). 


The results of the rotational derivative calculations depends upon the amount 
of detail used in the geometrical description of a shape. This effect is illus- 
trated in Figures 20 through 23. These figures contain the results of the 
ae calculations for a cone and for a wedge. The geometrical representa- 


tions used are shown in Figures 20and 21. The results obtained with different 
amounts of geometry detail in the longitudinal direction (AX/body length) are 

shown in Figures 22 and 23. These figures show that the program calculated 

values of Con are within less than one percent of the exact analytical 


q 
results when the shape is represented with longitudinal divisions of 10 per- 
cent of body length or less. It is evident from these results that the normal 
definition of a vehicle as used to give accurate static coefficients is also 
quite adequate for the dvnamic derivatives. 
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Figure 20. 20° Half Angle Cone. Graphical Representation of Various 
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The derivative Con _is much more difficult to calculate accurately since 
& 

it is caused by unsteady flow phenomena. Fortunately, this term is 

quite small for hypersonic vehicles in high-speed flight so approximate 

methods usually suffice. This stability derivative is the change in pitch- 

ing moment with rate of change of angle of attack and occurs at constant 

pitch angle. It is frequently referred to as the "plunge" derivative 

since it represents the effect of a change in vertical acceleration. 

The common definition of this derivative is 


ak 
vee” old) 


For a rigic vehicle this derivative arises from an aerodynamic time lag 
effect since the flow is not able to instantaneously adjust to changes in 
flight conditions. For low-speed airplanes the major component of this 
coefficient is caused by the fact that a change in the wing downwash 
field takes a finite length of time before it arrives at the tail. For 
tailless aircraft chere stillis a co _ Since the vehicle must accelerate 
a 
the air mass in its path as it accelerates. This parameter is also sub- 
ject to high~speed aeroelastic effects. 


No general theory is currently available for calculating CU _ The 
: a 
parameter itself is usually difficult to obtain experimentally since it 
cannot be separated from the total damping coefficient measured 
( Cn + C4 . References 60, 61, and ¢2 contain discussions of 
x q 
the most frequently used procedures for calculating Ca for wing- 


body-tail configurations. Reference 63 discusses the use of the 
rheoelectrie analogy for stability derivative determination. 


In the analysis for the damping derivative Ce. it is convenient to divide 
a 


the vehicle into separate components consisting of the bady, wing, and 
tail surface. The wing contribution may be found from the relationship 
below. 


(os) 


where 


1 
ome 
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Hi 


si exposed wing area 


S = reference area 
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K = interference factor for effect of wing in 


we presence of body (see Reference 64) 

Kaw = interference factor for effect of body in 
presence of wing 

a = mean aerodynamic chord of exposed wing 
c = reference chord for moment coefficients 
Gi = _ basic wing @ derivative as obtained from 

m. 

a Reference 60. 


The body contribution may be obtained by using the relatively simple 
results derived from slender~body theory. Reference 65 points out 
that although slender-body theory alone does not accurately predict 
the characteristics of nonslender configurations, the ratio of slender- 
body derivatives may be used to obtain reasonably accurate results 


as indicated below. 
Ca: 
Cc & 
My B a 
; slender body 


a 





@) 
3 
Sa 
93) 
ul 


body pitching moment derivative as calculated 
by the arbitrary body program for the actual 
shape involved 


a 
3 
© 
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center-of-gravity and area-centroid 
locations relative to the overall body 
length 
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S = bedy frontal aren 
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The contributions of a horizontal tail to Co. are given by the relationship 








@ 
‘ 2Q Sy -— == 
Cc = - s* ( ) (K +K “ty 2 (Se 
( me) gO Sua (ts Bt) (2 ACH 
where 
a (Ee); 
Q = tail effectiveness ratio = a (,). 
co G.)., 
Sp = tail exposed area 
= tail dihedral angle 
Pp 
(Cx ) = tail lift curve slope 
ry 
T 
King and Kar = tail lift interference factors 
Xep = tail length 
a (=2) = average upwash induced by the wing at the 
aZ\aV tail location (becomes negligible in the 


hypersonic range) 


The accuracy of the above methods for high angle of attack conditions 
or for very blunt lifting re-entcy bodies is not known, although Reference 
63 states that slender body theory is useful for such situations. 


The use of the above relatively simple approach for the calculation of 
the Ci derivatives would allow the determination of the 6 derivatives 

a 
in a similar manner (see Reference 65), 

Cc 
Y . 

C = (Cy } wen 

B slender body 
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: Volume 
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SECTION III 


PROGRAM ORGANIZATION 


The basic organization of the Mark III version of the Arbitrary-Body 
Program is shown in the diagram below. Each major part of the program 
could, with only minor modifications, operate as a completely separate 


computer program. 









EXECUTIVE PROGRAM 
(Ha8S} 









SLAB DELTA PROGRAM CARD PUNCH PROGRAM | 
{SLABD) (CARG) 





[AERODYNAMIC PROGRAM PGRAFHICS CONTROL 
(AERO, (GRAPIC) 
PICTURE DRAWING 
PROGRAM 
(PICTUR} 


The most important component of this system is the Aerodynamic Program 
(AERO). The general features of the logic flow for this part of the system 
are shown on Figure 24. Each subroutine is described in more detail in 
Volume I of this report (pages 16 through 25). « .: graphics parts of the 
system are described on pages 25 and 26 of Volume I. 








OUTPUT DATA 
PLOT'SR PROGRAM 
(PLOT) 
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start | 


He Ne ee ee ee ae see FRED See ero ome Suey ener 


AERO 


READ ELEMENT DATA 
TITLE CARD 







FRET! =] 


mest READ FORCE DATA [ag 
% TITLE CARD _ az 


[READ FLIGHT CONDITION CARD 
|M, H, FLAGS, NO. OF a-8, ETC.} 






READ FLIGHT ATTITUDE CARD | 
Qa, B, K, Se £ ETC. 


| GO TO FORCE PROGRAM 


SET NEW @ OR BCR 


FOR DERIVATIVES 


CHECK CONTROL | __ 
DEFLECTION —— 






CHECK FOR LAST a-8 


See A ST CURRED CURES GEES CEES GERM GERED GuETErEE <EICWES GuamnEStr GLUED GCrEEERD GEER GEESE GEEKS GREED Geewe GEE Gee? CRED GEE CED GONE: 


ge{ READ SKIN FRICTION CARDS] 


wel SET UP DERIVATIVE DATA 


- CALCULATE DERIVATIVES 


SURFACE DATA 


BBOREAD IN GEOMETR' 
j CALCULATE QUAL 













CONTROL OEFL 
(CONTRE 





CALCULATE FORCI 


ax-B DERIVAT! 
ELLIPTIC. INTEC 
= THRUST FORC 


Be) PRINT PAGE HEAL 












(SDATA} | PARAMETRIC CuBr ~] 
| ccf METRY DATA ag 
ets (ANALY 2) 


RILATERALS fagt>we 





ANALY 3 (DUMMY) 


SHOCK-EXPANSION 
(SHKEXP) 






CTION | 


$ FORCE) Fy 


aera 


CONTROL SURFACE 
FLOW_ SEPARATION 


| BLUNT-8ODY VISCOUS 
: (BLUNT) 

aaa ATMOSPHERE 

'eS (PLUNGE) | (ATMOS) 


RALS (ELP 1) | 
S (VECTOR) | 
INGS (HEADER) | 


<N 


[ ELUPTICAL SURFACE Le 
GENERATION (ANALY 1} 


FLOW SEPARATION (FLOSEP} LW. 





a SKIN FRICTION 
(SKINFR) 
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WRITE GEOMETRY DATA | 





- ON TAPE 8 
(PUNCH) 





a-| PRANDTL-MEVER EXPANSION 


(EXPAND) 


Ss 


TANGENT WEDGE 
(COMPR} 
'BLUNT-BODY NEWTONIAN + 
PRANDTL-MEYER (NEWTPM) 


SURFACE TEMPERATURE 
| (TEMP) | 


KEATING RATES 
(QC) 


Py 










FLUID PROPERTIES 
(ROMU) 






POLYt ‘“!ALS 
Pu.) 
THERMODYNAMIC DATA 
(DATA) 







Figure 24, AERO Program Flow Chart 


SECTION IV 


OPERATIONAL CONSIDERATIONS 2% * 


The Mark Ill version of the Hypersonic Arbitrary- ~Body Aerodynamic 
Computer Program System is written entirely in FORTRAN. Models. of 
the Mark Ill program are available for ase on the IBM 360, the IBM 7094, 
and the UNIVAC 1108 computers. For the Mafk Il program the IBM 360. 
model is considered to be the base-line program. The differences between 
the programs are of a minor nature and modifications necessary for 
operation on other similar computers should be easy to accomplish. The 
program also makes use of the Douglas version of the SC-4020 software 
system package to generate a plotting tape. This tape is then processed 
off-line by a Stromberg Carlson SC-4020 Data Recording System, The _ 
etaphics parts of the program may be converted by the user for operation 
with on-line graphics équipment such as the IBM 2250. 


OVERLAY STRUCTURE 


Because of the iarge size of this program it is necessary to use the overlay 
feature of FORTRAN. The overlay structure for the Mark III program is 
shown below. With this structure the program requires approximately 
105,.600 bytes of storage on the IBM 360 computer. The program will also 
operate or. an IBM 7094 (32K) with the requirement that eleven tape drives 
be available (including the standard input Tape 4, cutput Tape 6, and the 
SC-4020 output tape). 


os OVERLAY STRUCTURE 
MARK II! PROGRAM 


| MAIN ; 
TAPE 





PICTUR 
j_HEAOR 2 
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A 


" DECK. SET-UP AND. OPERATION - 
3 * te 


The overlay structure caitds for the IBM 360 computer are given below. 


ENTRY MAIN _ 
INSERT MAIN= — 
OVERLAY ALPHA 
INSERT AERO=, HEADER= 
OVERLAY BETA. 
INSERT SDATAz, ANALY l=, ANALY2=, ANALY3=, PUNCH= 
OVERLAY BETA 
INSERT CONTRL= 
OVERLAY PETA 
INSERT FORCES 
* INSERT IHCSERF 
* OVERLAY GAMMA 
_ INSERT SHKEXP=, COMPR=, EXPAND=, NEWTPMs, TEMP=, QC, 
.POLY=,ROMUs, CONES 
. sOVERLAY EPSILN 
Re INSERT FLOSEP= 
gt aye OVERLAY EPSILN. 
INSERT SKINFR= 
. OVERLAY EPSILN 
INSERT BLUNT= 
OVERGAY GAMMA 
INSERT ATMOS= 
OVERLAY BETA. 
"INSERT PLUNGE=, ELP 1s “- 
OVERLAY BETA — 
INSERT VEC TOR= 
OVERLAY ALPHA 
INSERT GRAPHIC= 
INSERT NXV=, XSCALV=, NHCAWAU, NHCAWAO, NHCAWAP * 
INSERT YMODV=, XMODV=,NHCAWBN * 
INSERT NHCAWFC, VXAXV=, APLOTV=, APRNTV=,BNBCDV=, * 
BRITEV=, DOTLNV= * 
INSERT ERMRKV=, ERRLNV=, ERRNLV=, GRIDIV=, LABLV=, LINEV=, * 
LINRV= * 
INSERT NOFRV=, NONLNV=, PLOTV=,SETCIV= * 
INSERT SCSETV=,ARGQ= * 
OVERLAY DELTA 
INSERT PICTUR=, HEADR2= 
OVERLAY DELTA 
INSERT PLOT= 
OVERLAY ALPHA 
| INSERT SLABD=, TTABLE =, QINT= 
| 


BEE Sia tent 


OVERLAY ALPHA 
INSERT CARD= 


Note: Those cards marked with an asterisk (*) contain the 
SC-4020 routines. The asterisk is not actually punched 
on the cards but is just used in this report to indicate 
the SC~-4020 insert cards. 
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oh BO 
Bahl 


$IOB 
$SETUP 
$RESTORE 
$* 


$EXECUTE 
$IBIOB 


HABS 


$ORIGIN 


AROA 
AROB 


$ORIGIN 


AROC 
AROD 
AROE 
AROF 
AROG 


$ORIGIN 
AROH 

$ORIGIN 
AROI 

$INCLUDE 


$ORIGIN 
| AROJ 
{ AROM 
: ARON 
j AROO 
i AROP 
AROR. 
AROS 
AROT 
AROU 
AROV 


The deck set-up for the IBM 7094 model of this program is outlined below.. . 


as required 
as required 
required on Douglas system 


as required to give on-line message to machine 
operator (any text in card columns 3~72} 


IBJOB 


Main program 


ALPHA 
Subroutine AERO 
HEADER 
BETA 
Subroutine SDATA 
ANALY! 
ANALY2 
ANALY3 
PUNCH 


BETA 

Subroutine CONTRL 
BETA 

Subroutine FORCE 
FERF 


GAMMA 


Subroutine SHKEXP 
COMPR 
EXPAND 
NEWTPM 
CONE 
TEMP 
QC 
POLY 
ROMU 
DATA 
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IBM 7094 Deck Set-Up (Continued) 


$ORIGIN 
AROL 
$ORIGIN 
AROK 
SORIGIN - 
AROQ 
$ORIGIN 
AROW 


$ORIGIN 
AROX 

- AROY 
$ORIGIN 
AROZ 


$ORIGIN 
GRPA 


$INCLUDE 


$ORIGIN 
GRPB 
GRPC 
$ORIGIN 
GRPD 


> $ORIGIN 


SLBA 
SLBB 
SLBC 


$ORIGIN 
CARD 


$DATA 


EPSILN 


Subrowvtine 
EPSILN 
Subroutine 


EPSILN 
Subroutine 


GAMMA 


Subroutine 


BETA 


Subroutine 


BETA 


Subroutine 


ALPHA 


Subroutine 


SKINFR 
FLOSEP 
BLUNT 
ATMOS 


PLUNGE 
ELP1 


VECTOR 


GRAPIG 


All SC-4020 subroutines required by program except 


those that must be in the main link. 


DELTA 


Sub routine 


DELTA 
Subroutine 
ALPHA 


Subroutine 


ALPHA 


Subroutine 


PICTUR 
HEADER2 


PLOT 


SLABD 
TTABLE 
QINT 


CARD 


All program input data cards 


Ind of File (7 and 8 punches in card column 1) 


$IBSYS 
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TAPE ASSIGNMENTS 


The operation of this program requires the availability of 11 logical units 
(besides the systems units). On the IBM 7094 these will be magnetic tape 
units, On the IBM 360 computer some or all of these units will be on disk 
storage. In either case, the logical units are referred to as "Tapes" 
throughout this discussion. The units required and their use in the program 
are listed below. 


Unit Mode Program Usage 


1 Binary Storage of aerodynamic coefficients 
for summation. 


3 Binary Storage of control fore-surface 
: geometry data. 


4 Binary Storage of excess of quadrilateral 
element data when number of 
elements is greater than ISIZE, 
Also used to store geometry data 
for control-surface flap. 


5 BCD Standard system input tape unit. 

6 BCD Standard system output tape unit. 

- Binary Standard system punch tape unit. 

8 BCD Storage of geornetry data in surface- 


element form. 


9 Binary Storage of aerodynamic data to be 
plotted by PLOT routine. 

10 ' Binary Storage of counter of aerodynamic 
data saved for plotting by PLOT 
routine. 

ll Binary Storage of control-surface geometry 


in the flap-deflected position. 


16 Binary SC-4020 output tape for the Douglas 
version of the SC-4020 system. 
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//G60.FTO3FOOL 
// 
//G0.FTO4FOO1 
// 
//GO.FTQ8F00L 
‘/ 
//GGeFTLIFOO:i 
Vi 
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Tape Assignments (Continued) 


For operation on the IBM 360 computer it is recommended that the following 
units be on magnetic tape: 3, 4, 8, and 11. These tapes may be defined 
as follows: 


DO DSNAME=8FTO3yUNITH=TAPEA sfISP=(NEW, DELETE), 
DC B=(RECFM=V, BLKSIZE=800,0EN=2, TRICH=C) 

DD DSNAME=EFT04, UNI T=TAPEB yy DISP=(NEWSDELETE) », 
OC B=(RECFM=V_BLKSIZE=800,DEN=2, TRICH=C ) 

DD DSNAME=TFTO8,UNIT=TAPEB ,DI SP=(NEWSKEEP) » 
OCB={ TRTICH=ET,RECFM=F, BLKSIZE=84,DEN=2) 

DD DSNAME=EFTL1,UNET=TAPESB DI SP=(NEW,DELETEL), 
DCK=(RECFM=V,BLKSIZE=800,UVEN=2 »TRICH=C) 


SC -4020 SYSTEM 


This program makes use of the Stremberg Carlson 5C-4020 Data Recording 


System. The SC-4020 software system used by this program was originally 
developed by North American Aviation, Inc. The Douglas version of this 


Sone en an mmienainensdndett tment himentiiahame ss nhbememaened 


$INCLUDE 
$INCLUDE 


ee ee 


$ INCLUDE 
$ INCLUDE 
$INCLUDE 
$ INCLUDE 


system contains a number of revisions. Some of these revisions were made 
to facilitate the use of the overlay features of FORTRAN, 


Because of the large core-storage requirements of the Hypersonic Arbitrary- 
Body System, it is important that all system subroutines be included in the 
proper overlay link. On the Douglas version of the S5C~-4020 system two 
subprograms must be overlay link 0. Other versions of the SC-4020 system 
may have eimilar restrictions. For the IBM 7094 model of the program 
these subroutines are PLOl and SCOUTQ. Those subroutines required in 
the overlay structure for the IBM 360 model of the program are indicated 
by the asterisk on page 116. All other SC-4020 routines on the IBM 360 
must be in the main link. The following cards may he used to include the 
appropriate SC-4020 subroutines in the proper part of the overlay structure 
for the IBM 7094 model of the program. 


APLOTV,APRNTV, BUTTV, BNBCDV, BRITEV,CO4020,CH4020 
DOTLNV,ER4020, ERMRKV, ERRLNV, ERRNLV, PR4020,GRIDIV 
HD4020,HOLDIV,HOLLV, 1D4020,INCRV, INTBCD,LABLV 
LE&O20,LINEV,LINRV,NOFRV,NONLNV, PL4020,PLOTV 
PRINTV,SCBSR,SETCIV,SETMIV,SMXYV,STOPTV, TF&620 
T14020,TPNUMV,WP4LO20,WR4020,XAXISV,XMODQ,XSCALQ 
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APPENDIX A a. oe 
PROGRAM LISTINGS AND FLOW CHARTS 


This part-of the reoort contains the listings and flow charts for the Mark Il 
Mod © version of the Hypersonic Arbitrary-Body Aerodynamic Camputer 
Program System: Also included.are the symbal lists defining the variables 


used in each subroutine, A combined list of all program variabies is included: 
in Appendix B, 


The symbol lists are divided into five fields which are described as follows: 


(i) The first field contains the symbol 


(ii} The second field contains the letters I, Ju, ox R, denotiag integer, 
logical, or real variable respéctively. 


(iii) The third field contains the letters A, C, D, ‘U, denoting argument. 
list, common, dimensioned, or undimensioned, variable respectively. 
The hierarchy of the above letters ig A, G, B, wv. 


{iv) The fourth field contains the definition of the symbol. 


(v) The fifth field contains the name of the subroutine in which-the 
symbol occurs. 


we ee BR mae oe 


cater et 8 tee 


INDEX TO SUBROUTINE LISTINGS AND FLOW CHARTS 


Routine 
AERO 
HEADER 
~ SDATA 
ANALY 
ANAY2. 
ANALY3 
PUNCH 
CONTRL 
FORCE 
SHKEEXP 
FLOSEP. 
SKINFR 
CQMPR 
EXPAND 
NEWTPM 
CONE 
BLUNT 
TEMP. 
QC 
POLY 
ROMU 
DATA 
ATMOS 
PLUNGE 
ELP1 
VECTOR 
GRAPIC 
PICTUR 
HEADR2 
PLOT 
SLABD 
TTABLE 
QINT 
CARD 


Name. - 


‘Main Program (Executive Program) 


_ Aerodynamic Program Coiitrol 


Pagé Header Stibroutine 


Surface ‘Data Subprogram 


Elliptical Cross-Section Subprogram 


Parametric Cubic Subprogram 
Analytical Shape Subprogram (Dummy) 
Blement Data Write Subroutine 
Control-Surfacé-Deflection Subprogram 
Force-Calculation Subprogram 
Shock-Expansion Subprogram 
Flow-Separation Subprogram 


Skin-Friction Subprogram 


‘Compression Subprogram 


Prandtl-Meyer Expansion Subprogram 
Blunt-Body Newtonian + Prandtl-Meyer 
Conical Flow Subprogram 
Blunt-Body Viscous Subprogram 
Temperature Subprogram 

Convective Heating Subprogram 
Polynominal Subprogram 

Equilibrium Air Properties Subprogram 
Block Data Subprogram 

Atmosphere Subprogram 
Plunge-Derivative Subprogram 
Elliptical-Integral Subprogram 
Thrust- Vector Subprogram 

Graphic Executive Subprogram 
Picture-Drawing Program 

Graphic Header Subprogram 

Output Data Plotter Program 

Slab Delta Program 

Triple Interpolation Subprogram 
Quadratic Interpolation Subprogram 
Card Punching Subprogram 


HABS  A-3. 
AROA = A-IL 
AROB A-47 
ARQOG  A-51 
AROD  A-75 
ARQE A-93 
AROF A-113 
AROG = A-l17 
AROH  A-123 
AROL A+133 
AROSJ As175 
AROK  A-187 
AROL A-217 
AROM A-247 
ARON A-257 
AROO } A-267 
AROP A-279 
AROQ = A-~285 
AROR A-293 
AROS A-309 
AROT A-319 
AROU A-323 
AROV A-333 
AROW } A-337 
AROX A-345 
AROY A-363 
AROZ A-369: 
GRPA A-377 
GRPB A-381 
GRPC A-413 
GRPD  A-4l7 
SLBA A-437 
SLBB A-455 
SLBC A-463 
CARD A-467 


‘Se 1, 





MAIN PROGRAM (DECK HABS) 


as 


c. 


Algorithm Peg 


This. routine, the main routine for the Mark III version. of the 
Hynersonic Arbitrary-Body Aérodynamic Computer Program 
System, acts as the system monitor. It controls and initiates 
éxecution. of the five principal components of this. systern. These 
are the Aerodynamic Force Analysis Program (Program AERO 
and its associated subroutines), the Graphics Prograrn Options 
(the Picturé Drawing Program. and: the Output Data Plotter Pro- 
gram), the Auxiliary Geometry Generation Program (Slab Delta 
Program), andthe card punch routine. 


The first action of the progtam. is for this: Executive routine to 
read in «.° System Control Card (Type 0). This card controls 
the selection of each program option and the order in which the 
options are to be used. Up to twenty different program phases 
may be used, and any given option may be used several times, 
If an error has been detected in any of the program subroutines 
control will be returned: to the Executive Program with the 
ERROR flag set to the appropriate value, Ifthe error is of the 
non~fatal type the Executive Main Program will start reading 
in data cards until it finds ‘a card with a ‘Type number of 99. 
The pr’ Jram will then attempt to execute the next phase option. 
Input/Output 


System Control Card (Type 0) 
Program Flow monitoring information 


Error 

An error condition occurs when an error is returned by a sub- 
routine, or when the System Control Card contains az illegal 
value. 

Subroutines Required 

AERO, GRAPIC, SLABD, CARD 

Argument List 

Not applicable 

Length 


2170 bytes 
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2. SUBROUTINE AERO (DECK AROA). 


This routine controls the flow of calculations within the Aerodynamic 
Force Caiculation Option. 


a. 


c. 


ox ner eater the RPT cs wih tn tt SOT OO OEE LCE CC: a 


Algorithm 


The first operation of this routine is ta rewind the necessary 
program tape units. it then reads the Element Data Title Card 
(Type 1) and calls the Surface Data routine for the geometry 
calculations. Upon return from. Surface Data (SDAA): the Force 
Data Title Card (Type 8) is read and if IRET1 is zero the Force 
Data Title Card will be read next; This is then followed by the 
Flight Condition Card, the Center of Gravity Data Card, and if 
required, the Skin Friction Data Cards and Coefficient Incremenit 
Data Cards, The program then starts the cycle to calculate 
force data at each flight attitude. This is initiated by the reading 
of the Flight Attitude Data Card: The necessary. data arrays 

are initialized, derivative data established, and the Force cal~ 
culation routine called. 


Upon return from the FORCE routine the derivative cycle data 
is set and control returned to the FORCE routine if required. 
If the control surface is to be deflected the CONTROL routine 
will be called before the FORCE routine is called, If required 
the Plunge Derivative and Thrust Vector routines wit also be 
catled, 


Input/Output 


Element Data Title Card (Type 1), Force Data Title Card (Type 8}, 
Flight Gondition Card (Type 9), Center of Gravity Data Card 

(Type 10), Skin Friction Data Cards (Type 11), Flight Attitude 
Data Cards (Type 12), and Coefficient Increment Data Cards 

(Type 13). Output data consists of the vehicle summary force 
data, 


Error 


An error condition occurs when input card Type number is in error. 


Subroutines Required 

SDATA, CONTRL, FOAGE, PLUNGE, VECTOR, HEADER 
Argument List 

None 

Length 

15396 bytes 
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%* Bs 


c. 


e. 


a RL 


Algorithm 


This routine provides the title at the top-of each page of the cutput 
and advances the page counter, 


Input/Output 


Program header is printed at top of page on output Tape 6, 


Error 
None 
Subroutines Required 


None 


Argument List 


None 
Length 
342 bytes 
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SUBROUTINE SDATA (DECK AROC) 


This subroutine prepares. geometry data for use by the rest of the 
program. 


a. 


Cc. 


e. 


Algorithm 

The Element Data Control Card (Type 2) is read and control passed 
to one of the other geometry routines (ANALY1, ANALY2, ANALY3) 
if required. The surface element data is then read (either from 
input Tape 5 or from the geometry storage Tape 8) and converted 
to quadrilaterals, These data are stored in core for the first 

300 elements and on Tape 4 thereafter. All geometry data for 
control surface components are stored on Tapes 3 and 4. 
Input/Output 

Element Data Control Card (Type 2), Elernent Data Input Cards 
(Type 3). 

When PRINTS is equal to 1 the Input Surface Element Data along 
with the direction cosines and centroid coordinates of each 
quadrilateral element are printed on output Tape 6, 

Error 

An error condition occurs when an input card type number is wrong. 
Subroutines Required 

ANALYI, ANALY2, ANALY3, HEADER 

Argument List 

(PRINTS, SYMFCT, ISIZ, IORIEN, IGTYPE) 

Length 


13264 bytes 
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5, SUBROUTINE ANALY] (DECK AROD) 


This subprogram prepares surface-element data points for circuit. »~ 
elliptical cross-sections. 


a. 


ee, 


Algorithm 

The Ellipse Generation Control Card (Type 4) is read as are all 
the Cross-Section Data Cards (Type 5), The input cross-section 
data are stored in the appropriate arrays. The program then 
takes each cross-section and calculates the surface element data 
as directed by the input parameters. This information is stored 
on Tape 1] and at the appropriate time is transmitted to the 
PUNCH routine, recording onto the geometry storage tape. 
Input/Output 

Ellipse Generation Control Card (Type 4), and Cross-Section 
Data Cards (Type 5). 

No output data are produced. 

Error 

An error condition occurs when input card type numbers are wrong. 
Subroutines Required 
PUNCH 

Argument List 

None 


Length 


8552 bytes 
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6. SUBROUTINE ANALY2 (DECK AROE) 


This routine prepares surface element data from parametric cubic 
boundary patch data. 


a. 


Algorithm 


The routine first reads ir. the Parametric Cubic Title Card 
{Type 6). The boundary curve data for a patch is then read in 
(Type 7 cards). A maximum of 20 points per boundary are per- 
mitted, The boundary curve coordinate arrays are set up and 
the arc length along each boundary calculated. The tangent 
vectors at the corner points and the related end-point derivatives 
are calculated. The constants for the boundary curve equations 
are then determined. The routine then starts a cycle to gener- 
ate element data points on the surface of the patch. This requires 
the determination of the blending functions, the related points on 
the boundary curves, and finally, the solution of the surface 
equation. The element data are transmitted to the PUNCH 
routine where they are recorded onto Tape 8. 


Input /Output 

Parametric Cubic Title Card (Type 6), Parametric Cubic 
Boundary Data Cards (Type 7). 

No output data are produced. 

Error 

An error condition occurs when input card type numbers are wrong. 
Subroutines Required 

PUNCH 

Argument List 

None 


Length 


6688 bytes 
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SUBROUTINE ANALY3 (DECK AROF) 
This is a dummy subroutine provided for future use. This routine 
may be replaced with a routine similar to either the Ellipse Gener- 
ation routine or the Parametric Cubic routine. 
a. Algorithm 
b. Input/Output 
c. Error 
d. Subroutines Required 
e. Argument List 
None 


f. Length 


124 bytes 
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8 SUBROUTINE PUNCH (DECK AROG) 

This routine writes element data card images on Tape &. 

a. Algorithm 
If IPRINT is equal to 1, each card image will be written on 
output Tape 6. The element data are recorded on Tape 8 
in exactly the same form as normal input surface elernent 
data (Type 3). Each card is given a sequence number and 
the number of records written on Tape 8 is furnished to the 
calling routines. 


b. Input/Output 


Element data are recorded on Tape 8 and also on the standard 
output tape if required. 


ec. Error 
None 

d. Subroutines Required 
None 

e. Argument List 


| (X1. Yl, 41, NSTATI, X2; Y2, 22, NSTATS,. SECT, TYPE, 
LINE, SEQ, LAST, IPRINT, NREC) 
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SUBROUTINE TONTRL (DECK AROH) 


This routine changes the geometry data for control surfaces to the 
proper deflected position. 


a. 


e. 


Algorithm 

The hinge-line coordinate data are determined from the first 
two items in the NX2 data array. The required transformation 
angles are then calculated. The element data for the control 
flap in the undeflected position are read from Tape 4. The 
rotation matrix and the control deflection matrix are then 
applied to the geometry data. These data are then rotated 
back to the original hinge-line centered coordinate system, 
Hinge moment factor data are calculated and stored in the 
XCENT2 data array. The new geometry data for the deflected 
flap are stored on Tape 11. 

Input/Output 

None 

Error 

An error condition occurs when the total number of elements 
on the control surface is greater than ISIZ (300 for the Mark II 
program). 

Subroutines Required 

None 

Argument List 

(DELTAE, ISIZ) 

Length 


2376 bytes 
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10. SUBROUTINE FORCE (DECK AROI) 


This routine determines the pressure coefficients on each quadrilateral 
element, resolves the force in the required body axis system, and sums 
the contributions of each element to give the vehicle's six aerodynamic 
coefficients. 


a. 


Cc. 


Algorithm 


First the necessary starting constants and conditions are set up. 
The free-stream properties are then determined either frem the 
Atmosphere subprogram or frorn the ideal gas equation of NASA 
TR 1135. The loop for calculating and summing up the forces on 
each element is then started. The direction cosines of the velocity 
vector are caiculated and the quadrilateral element data are read 
from core or from tape. Symmetry requirements are then checked 
and the signs changed if required. The velocity components with 
vehicle rotation and the impact angle are calculated. If the impact 
angle is less than zero then the element is in a shadow region; if 


not, it is in an impact region. The proper force calculation method is 


then used to determine the pressure coefficient, The force sum- 
mation method to meet symmetry requirements is selected, the 
six aerodynamic coefficients calculated, and summation of the 
coefficients is accomplished. If required, the Skin Friction 
Subprogram is called for the viscous calculations. 


Input/Output 


If PRINT is equalto 1 the detailed force contributions of each 
element will be printed. 


Error 


An error will occur when the arccosine of the angle »etween the 
velocity vector and the unit normal is greater than I, 


Subroutines Required 


ATMOS, NEWTPM, COMPR, SHKEXP, EXPAND, FLOSEP, 
HEADER, SKINFR, BLUNT, CONE 


Argument List 


(ALPHA, BETA, CPSTAG, SREF, SYMFCT, XCG, YCG, ZCG, 
MACH, SPAN, MAC, J, QRP, ALT, PRINT, CAI, CNI, CYI, 
CLLI, CLMi, CLNI, ETAC, NS, IMPACT, IPRINT, IFIRST, 
PSTAG, TSTAG, RENO, ISIZ, ENPM, QQINF, ISHAD, IORIEN, 
IMPACI, JSHADI, IDERIV, V, IGTYPE, DELTAE, SWEEP, 
RETRAN, HML, HMR, PFS) 


Length 
13156 bytes 
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ll. SUBROUTINE SHKEXP (DECK AROJ} 


This subroutine performs a shock expansion analysis along a stream- 
wise strig of elements. 


a. 


ony 


&. 


lzorithm 
The element is first checked to see if it is the first element in a 
strip. Ifitis, the local properties on it are calculated by the 
appropriate methed and saved for use on the next element, Fer 
the next elernent in the strip the turning angle is calculated and 
either the compression or the expansion routine is used to deter. 
mine its local properties (pressure coefficient, Mack sid 
and temperature}. : 
Input/Output 
Nene 
Error 
An error condition occurs when the number of initial strip 
elements is greater than 100, and when the wrong initial element 
method has been input. 
Subroutines Required 
COMPR, EXPAND, CONE 
Argument List 
(ORIEN, N, M, IPRCK, NX, NY, NZ, DELTA, PFS, MACH, 
IiGTS, DELTAR, IMPACI, CPSTAG, TFS, CP, G, ISHADI, 
IFIRST, LL, IPRINT, IGTYPE, RHOFS, AFS, VIS, V, RENO, 
ISMODE) 
Length 


4472 bytes 
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12. SUBROUTINE FLOSEP (DECK AROK) 


This subroutine has the task of determining the effect of flow separation 
caused by the deflection of a control surface. 


a. 


c. 


Algorithm 

The flow separation calculations are performed in three cycles, 
The first cycie is with the control defiected and determines if 
seperation will occur. The second cycle determines where the 
Separation occurs and the pressure changes due to separation. 
The third cycie caiculated the finel pressure, using the normal 
program input methed in combination with the flow separation 
increments. 

Input /Gutput 

None 

Error 


An error condition occurs when the number of streamwise elements 
for the control fore-surface and flap is greates than 125. 


Subroutines Required 

TEMP, COMPR, SHKEXP 

Argument List 

(IGT, LL, N, M, NX, NY, NZ, XCENT, YCENT, ZCENT, AREA, 
XPA, YPA, ZPA, IGTYPE, IORIEN, IPRCK, DELTA, PFS, 
MACH, DELTAR, IMPACI, CPSTAG, TFS, CP, G, ISHADI, 
IFIRST, ISBP, IFSCY, IGTS, IMPS, ISHS, L, XLE, DELTAE, 
SWEEP, CPNIN, DELCPC, KATACH, XSEPP, ISPNT, IMPACT, 
ISHAD, RETRAN) 

Length 


19208 bytes 
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13, SUBROUTINE SKINFR (DECK AROL) 


This routine calculates viscous forces for both laminar and 
turbulent flows including viscous-interaction and planform 
effects. Also, the viscous induced pressures are determined, 


a. Algorithm 


The basic program constants and control flags are first 

established. The routine then starts a major DO-loop to 

calculate the viscous forces on each of the skin-friction 

surfaces, This involves the following steps: Caiculate 

the surface planform geometry, local flow conditions, ~ 
surface temperature, viscous~interaction effects, laminar : 
and turbulent visc-us forces, and summation of forces for 

either laminar or turbulent as specified. If desired, the 

induced pressure increment due to viscous -interaction is 

also calculated. At the user's option, the skin friction may 

be calculated by the reference temperature, reference 

enthalpy, or theSpalding-Chi methods. 


b. Input/Output 


If IS(I, 7) = 1, skin friction data for each surface will be 
printed. 


c. Error 


An error condition will occur if the input Mach number is 
subsonic. 


d. Subroutines Required 
TEMP, CONE, COMPR, EXPAND, NEWTPM, HEADER 
e. Argument List 


(ALPHA, CA, SREF, SKIN, NS, ALT, MACH, CPSTAG, TFS, 
PFS, AFS, RHOFS, IFIRST, VIS, ISIZ, CN, IGTYPE, DELTA) 


f. Length 


10,030 bytes 
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14, 


SUBROUTINE COMPR (DECK AROM) 


This routine calculates the local flow properties using conventional 
oblique-shock relationships (NACA TR 1135). 


a. 


Algorithm 

First the constants in the oblique equation are calculated. A 
check is made for shock detachment, and if the shock is not 
detached the three real roots for the cubic are found. If the 
shock has detached, the conditions will be calculated by the 
NEWTPM routine. The proper root is selected and local flow 
conditions are caiculated. - 

Input/Output 

None 


Error 


An error condition will occur if a negative value is found for 
sine theta squared, Set it to 0 and the program will continue. 


Subroutines Required 
NEWTPM 
Argument List 


(ANGLE, MER, IPRINT, CPSTAG, TFS, PFS, ISDET, IFIRST, 
CP) 


Length 


2794 bytes 
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SUBROUTINE EXPAND (DECK ARON} 


This routine calculates the local flow conditione on a surface by using 
conventional Prandtl-Meyer relationships (NACA TR 1135}. 


a. 


b, 


Algorithm 

Calculate the Prandtl-Meyer angle after the expansion. Check 
if the flow is to be cornpressed to subscnic conditions or if it 
has expanded to an infinite Mach number (100). An iterative 
procedure is used to deterrnine the flow conditions after the 
expansion fan. 

Input/Output 

None 

Error 

None 

Subroutines Required 

None 

Argument List 

(ANGLE, MER, IPRINT, ISDET, CP) 

Length 
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SUBROUTINE NEWTPM (DECK AROO} 


This subroutine calculates the pressure coefficients on a surface by 
the biunt-body Newtonian + Prandtl-Meyer method. 


a, 


Algorithm 

This first section of this routine performs an iteration to find 
the matching point Mach nurmber. A Prandtl-Meyer expansion is 
then calculated from the matching point condition to the locai 
element slope using the EXPAND routine. Finally, the pressure 
coefficient and the locz! flow properties are calculated. 
Input/Output 

None 

Error 

None 

Subroutines Required 

EXPAND 

Argument List 


{ANGLE, EMN, CP, ETAC, IPRINT, MER, CPSTAG, TFS, 
P¥FS, ISE, IFIRST) 


Length 


2862 bytes 
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17, SUBROUTINE CONE {DECK ARCP} 


ae 


a. 


This subroutine sclves fer the local properties about 
a cone in supersonic flow using empirically derived 
equations. 

Algorithm 

Calculates the shock normal Mach number, surface 
pressure coefficient and, at theuser's optiou the 
following local flow properties on the cane surface: 
pressure, density, temperature, velocity, speed af 
sound, Mach number, viscosity, and Reynolds num- ._ 
ber per foot. Solutions are empirically derived for 
a calorically perfect gas with ratio of specific heats 
equal to 1,40. 

Input/Output 

None 

Error 

None 

Subroutines Required 

None 

Argument List 

(ANGLE, CP, ISDT) 

Length 


1256 bytes 
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18, 


ie 


SUBROUTINE BLUNT (DECK ARCQ) es kez 


This routine calculates the viscous forces on a blunt 
_body including low density effects. 


a. Algorithm 
Checks for ideal or real gas option, then calculates 
local properties behind normal shock. Determines 
local viscous forces and calculates low density 
viscous~interaction effects. 

b. Input/ Cutput 
IPRINT = 1, pertinent local and free-stream properties 
and viscous force coefficients will be printed. (This is 
intended for checkout only and IPRINT must be set 
within the program. ) 

ce, Error 
None 

d. Subroutines Required 


None 


e. Argument List 


(PFS, MACH, TFS, VIS, RHOFS, RB, RENO, TAU, IVISIN) 


f. Length 
1596 bytes 
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19. SUBROWTINE TEMP (DECK AROR) 


This routine uses an iterative procedure to calculate the surface 
equilibrium temperature for either an ideal gas or a real gas. 


a. Algorithm 


Calculate local and recovery enthalpy and general constants. 
Check for type of calculation desired {ideal or real gas, tem~- 
perature input or calculate), proceed with iteration and de- 
termine local convective heating rate, reference Reynolds 
number and compressibility factors. Print local skin friction 
data and temperature iterations if required. - 


b. Input/Output 


IPRINT = 1, temperature iterations and local skin friction 
data will be printed. 


IPRINT = 2, only local skin friction data will be printed. 


c, Error 


a None 


d, ..Subroutines Required 

QC 
e. Argument List 

(EL, TR, RE, TS, NW, MER, IPRINT, RT) 
f. Length 


3620 bytes 
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20. FUNCTION QC (DECK AROS) 


This routine calculates the aerodynamic heating at the given wall 
temperature. 


a. 


e. 


Algorithm 

Tests for laminar or turbulent flow, for reference method 
or Spalding-Chi, and for ideal or real gas. Calculates 
convective heating rate and sets certain quantities in 
common. 

Input/Output 

None 

Error 

None 

Subroutines Required 

ROMU (three entries, ROMU, ROW, ENTHAL) 
Argument List 

(TW) 

Length 
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21, FUNCTION POLY (DECK AROT) 


This routine generates an N-th order polynomial, 


a. 


Cs 


Algorithm 


Polynomial is evaluated for the input order and coefficient 
array at a specified starting value. 


Input/Output 

None 

Error 

None 

Subroutines Required 
None 

Argument List 

(A, I, HX, N) 
Length 

444 bytes 
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22. FUNCTION ROMU {DECK AROU) 


. This routine calculates various equilibrium air real gas 
properties. Has three entries; ROMU, ENTHAL, and ROW. 


a. Algorithm 


ROMU determines the density-viscosity product as a 
function of the input enthalpy and pressure. 


ENTHAL determines the enthalpy as a function of the 
input temperature and pressure. 


ROW calculates the density as a function of the input 
enthalpy and pressure. 


b. Input/Output 
None 
c. Error 
None 
d. Subroutines Recuired 
POLY 
e. Argument List 
ROMU (HS, P2), ENTHAL (TW, PW), ROW(HS, P2) 
f, Length 
2478 bytes 
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23. BLOCK DATA (DECK AROV) 


This routine initializes data arrays into common required in 
calculating equilibrium air real gas properties. 


a. 


Algorithm 


Data arrays are initialized at time of compilation into 
labelled common /PROP/. 


Input/ Output 

None 

Error 

None 

Subroutines Required 
None 

Argument List 

None 

Length 

1080 bytes 
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24. SUBROUTINE ATMOS (DECK AROW} 


This routine caiculates the atmospheric properties using the 1962 
U. 5. Atmospheres, 


2. 








| 
| 


Algorithm 


Set up arrays and constant values. Calculate atmospheric 
properties assuming én invérse square gravitatienal field, 


Input/Catput 

Error 

None 

Subroutines Required 
None 
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25. SUBROUTINE PLUNGE (DECK AROX} 


This routine calculates Cm; and CY, for wings, bodies, and tails, 


Be 


Cc. 
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Algorithm 


Read the Plunge Derivative Control Card (Type 14). Check 
correctness of input data flags. Read in plunge derivative 
data (Card Types 15, 16, 17, 18, 19, 20, 21) as required, 
Select proper calculation method, determine KBW if required, 
and calculate final derivatives. 
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Input/Output 
Card Types 14, 15, 16, 17, 18, 19, 20, and 21 
Error 


An error condition occurs if the input flags are wrong, or if the 
card types are in error. 
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26. 


SUBROUTINE ELP1 (DECK AROY } 


Thie routine approximates the values of he elliptical integrals of 
the first and second kinds. 


2, Algorithm 


The approximation to a value of the elliptical integral of the 
first kind is given by 


4 ; 4 
K(k) = (a, tat... sam )t(b + byt... + bq) In 


where n= 1-k° 


aye 1, 386294361 bo = 0.5 

a, = 0. 0966634426 by = 0. 124985936 
a, = 0. 0359009238 b, = 0. 0688024857 
az = 0. 0374256371 b, = 0, 0332835534 
a, = 0.0145119622 b, = 0. 0044178701 


The approximation to the value of the elliptical integral of the 
second kind is given by 


z 4 4 1 
E(k) = (aj tant... tag )+(b > + byt... + ban’) In > 


where yn = 1 - a 
a =1,.0 b =0,.0 
° ° 
a, = 0. 4432514146 hy = 0. 2499836831 
a, = 0. 0626064122 by = 0. 0920018004 
a, = 0. 0475738355 b = 0. 0406969753 
a, = 0. 0173650645 by = 0. 0052644964 


b. input/Ontput 
None 
c. Error 
None 
d. Subroutines Required 
None 
e. Argument List 
{AK, AKK, E, NERROR) 
f. Length 
828 bytes 
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SUBROUTINE VECTOR (DECK AROZ) 


This routine converts input thrust vector data to coefficients and 
adds the results to the vehicle aerodynamic coefficients. 


a. 


Cc. 


e. 


Algorithm 

Set up initial conditions and constants. Read in a force 
vector and convert to force coefficients, Print the results 
if required, Continue to read in iorce vector data until 
LAST equals one. 

Input/Output 


Thrust Vector Data Cards (Type 22) 
Thrust Vector coefficient contributions are printed if IPRINT = 1. 


Error 

An error condition occurs if the card type number is wrong. 
Subroutines Required 

HEADER 

Argument List . 


(MACH, PFS, SREF, XCG, YCG, ZCG, SPAN, MAC, ALPHA, 
CD, CL, CA, CY, CN, BETA, LOD, CLM, CLL, CLN) 


Length 


2304 bytes 
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SUBROUTINE GRAPIC (DECK GRPA) 


This is the Executive routine for the graphics part of the program, 


a. 


Ce 


e. 


Algorithm 


Print that GRAPHIC OPTION HAS CONTROL and select the 
proper graphic routine depending upon the value of IPROG, 


Input/Output 

None 

Error 

None 

Subroutines Required 
PICTUR, PLOT 
Argument List 
(IPROG) 

Length 

536 bytes 
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SUBROUTINE PICTUR (DECK GRPB) 


This routine prepares an output tape for procession on the SC-4020, 
The result will be pictures of the vehicle with the selected viewing 
angles. 


a. Algorithm 


Read the Picture Drawing Program Element Data Title Card 
(Type 33) and the Element Data Control Card (Type 32). Read 
the surface element data either from Tape 5 or from Tape 8 

as directed. Read plotting instruction data (Card Types 34, 35, 
36, and 37). Set up starting constants for pictures. Read 
element data from Tape 3 using the same techniques as for 
SDATA and convert to quadrilaterals, Generate scale grids if 
required, Plot points and draw lines between the points as 
directed by the input data. Print the detailed element char- 
acteristics if PRINTS is equal to l. 


b. Input/Output 
Element Data Title Card (Type 31), Element Data Control Card 
(Type 32), Element Data Cards from Tape 5 or Tape 8 (Type 3), 
Picture Control Data Card (Type 34), Grid Data Card (Type 35), 
Scale Label Card (Type 36), and Plot Title Card(s) (Type 37). 
If PRINTS is equal to 1 the detailed element characteristics will 
be printed on Tape 6 (just as is the case for SDATA). 

c. Error 
An error condition occurs if a card type number is wrong. 

d. Subroutines Required 
HEADR2, SC-4020 routines 

e. Argument List 
None 


f. Length 


16180 bytes 
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SUBROUTINE HEADR2 (DECK GRPC} 


ae 


Cs 


f, 


Algorithm 


This routine provides the title at the top of each page of the 
output and advances the page counter. 
similar to the HEADER routine. 


Input/Output 


Program header is printed at top of page on output Tape 6. 


Error 


Nane 


Subroutines Required _ 


None 

Argument List 
None 

Length = 
350 bytes 
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SUBROUTINE PLOT (DECK GRPD) 


This routine is used to produce graphically plotted data as obtained 
frorn the aerodynamics part of the program. 


ae 


e. 


Algorithm 

Read in plotter control cards. As directed, read aerodynamic 
data from Tape 9. Prepare plot scales and grids. Plot data 

and connect data points as directed. 

Input/Output 

Data Source Control Card (Type 41), Vertical-Title Card (Type 44), 
Horizontal-Title Card (Type 45), Plotting-Grid Data Card (Type 45}, 
Plot Control Array Card (Type 47), and Horizontal-Lable Card(s) 
Type 48). 

Output plots are on the SC-4020 tape. 

Error 

An error condition occurs when the card type number is wrong. 
Subroutines Required 

None 

Argument List 

None 


Length 


10772 bytes 
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26. SUBROUTINE ELPI (DECK ARGY } 


This routine épproximates the values of the elliptical integrals of 
the first and second kinds. 


a. Algorithm 


The approximation to a value of the elliptical integral of the 
first kind is given by 


“3 4 4 1 
K(k) = (a, tamt... + an + (bo + bint... + bn ) Ins 


where 7 =1 - Ke 


ay 1, 386294361 by = 0.5 

ay > 0. 0966634426 by = 0, 124985936 
a5 = 0. 0359009238 by = 0. 0688024857 
a = 0. 6374256371 b3 = 0, 0332835534 
ay = 0. 0145119621 b, = 0. 0044178701 


The approximation to the value of the elliptical integral of the 
second kind is given by 


ae: 4 4 1 
E(k} = f@, tant. + a,n Jt(b tht... thn) ce 
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where nH = 1 - i” 
a = 1,0 b =0.0 
° Q 
a, = 0. 4432514146 b) = 0, 2499836831 
as = 0. 0626060122 b5 = 0,0920018004 
az 0, 0475738355 b3 = 0. 0406969753 
’ a, = 0.0173650645 - Dy = 0, 0052644564 
b&b. Input/Output 
E None 
c, Error 
[ None 
f d, Subroutines Required 
None 
e, Argument List 
4 {AK, AKK, E, NERROR) 
: { 7, Lereth 
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zi, SUBROUTINE VECTOR (DECK AROZ) 





This routine cenverts input thrust vector data to coefficients and 
adés the resulis to the vehicle acradynamic coefficients. 


a. Algorithm 

Set up initial conditions and constants, Read in a force 
vecter and convert to force coefficients. rint the results 
if required. Continve to read in force vecto. data until 


LAST equals one. 


b. Input/Output 


ere 


ants 


Thrust Vector Data Cards (Type 22) 
Thrust Vector coefficient contributions are printed if IPRINT = 1, 


c. Error 

An error condition occurs if the card type number is wrung. 
d, Subroutines Requirei 

HEADER 

Argument List , 


(MAGH, PFS, SREF, XCG, YCG, ZCG, SPAN, MAC, ALPHA, 
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28. SUBROUTINE GRAPIC {DECK GRPA) 
Uiis is the Executive routine for the graphica part of the program, 
a, Argorithm 
Print that GRAPHIC OPTION HAS CONTROL and select the 
proper graphic routine deperding upon the value of IPROG, 
b, Input/Outpur 
None 
é ce. Error 
None 
d. Subroutines Required 
PICTUR, PLOT 


e. Argument List 
(IPROG) 

f. Length 
536 bytes 
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%. SUBROUTINE PICTUR (DaCK GRPB) 


Te: 


This routine prepares an cutput tape tor procession on the SC-4020, 
The result «ili be pictures ot the vehicle with the selected viewing 
angles. 


a. Algorithm 


Read the Picture Drawing Program Element Data Title Card 
{Type 31) and the Element Data Control Card {Type 32). Read 
the sucface clement data either frora Tape 5 o» from Tape 8 

as directed. Read plotting instruction data (Card Types 34, 35, 
36, and 37). Set up starting constants fur pictures. Read 
element data from Tape 3 using the same techniques zs for 
SDATA and convert to quadrilaterais, Geuerate scale grits if 
required. Plot points and draw lines between the pcints as 
Gairected by the input data, Print the detailed element char-~- 
acteristics if PRINTS is equal to l. 


b. Input/Output 


Element Data Title Card (Type 31), Element Data Control Card 
{Type 32), Element Data Cards frorn cape 5 or Tape 8 (Type 3}, 
Picture Control Data Card (Type 34), Grid Data Card (Type 35), 
Scale Label Card (Type 34), and Plot Tithe Card(s) (Type 37}. 

If PRINTS is equai to | the detailed eleraent characteristics will 
be printed on Tape 4 (just as is the case for SDATA), 


An error condition scecurs if a card type number is wrong. 
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30, SUBROUTINE HEADR2 (DECK GRPC} 
a.  Algorivhm 
This routine provides the title at the top of each page of the 
output and advances the page counter. This routine is very 
similar to the HEALER routine. 
b. Input/Output 
Program header is prirted at top of page on output Tape 6. 
c. Error 
Nene 
dad. Subroutines Requir .d 
None 
e. Argument List 
Notie 
f. Length 
350 bytes 
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33. SUBROUTINE PLOT (DECK GRPD) 


This routine is used to produce graphically plotted data as sbtained 
from the aerodynamics part of the program. 


a. Algorithm 
Read in plotter control cards. As directed, read aerodynamic 
data from Tape 9. Prepare pilot scales and grids. Plot data 


and connect data points as directed, 


dD. Input/Output 


ewan perenne 


Data Source Contro} Card (Type 41), Vertical-Title Card (Type 44), 
Horizontal-Title Card (Type 45), Plotting-Grid Data Card (Type 45}, 
Plot Control Array Card (Type 47), and Horizontal-Lable Card(s) 
Type 48). 

Ourput plots are on the SC-4020 tape. 


Error 


" 
e 


An error condition occurs when the card type number is wrong. 
d, Subroutines Required 

None 
e, Argument List 

None 
f. Length 


10772 bytes ee 
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32. SUBROUTINE SLABD (DECK SLBA) 


This routine generates the element data for a simple slab delta vehicle. 


a. 


Algorithm 


Read in Slab Delta Title Card (Typ=s 50), and the Configuration 
Control Card (Type 51). Hf required read in the Thickness 
Correction Cards {Type 52 and 53}, Read ina Cross-Section 
Deta Card (Type 54) and calculate the element data for this 
X-station. Write the card images on the -cgular output Tape 6 
if required and also on the geometry storage Taye 4. 

Continus until all the X-station cards have been read. 


Input/Output 

Slab Delta Title Card (Type 50}, Slab Delta Sweep Card (Type 51), 
and the Slab Delta Station Data Card (Type 54). [If ITOC is 
equal to 1 the Thickness Correction Cards (Type 52 and 53) 
are also inpvt. 

The card images cf the element data are written on the normal 
output Tape 6 if IPR{NT equals 1, and also on the geometry 
storage Tape 8. 

Error 

An error condition occurs when 2 card type number is wrong. 
Subroutines Kequired 

TTABLE 

Argument List 

None 


Length 


7584 tytes 
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33. 


SUBROUTINE TTABLE (DECK SLBB) 


This routine performs the interpolation to fird tne thickness correction 
factors for the Siab Dei:a Routine. 


a. 


Cc. 


e. 


Aigorithia 

Search fcr the proper points in the data table to be used in the 
interpolation. Call on the quadratic interpolation routine, 
QINT, to obtain the interpolated value, 

Input /Outpet | 

None 

Error 

None 

Subroutines Required 

QINT 

Argument List 

(A, B, C, D, R, G, Gl) 

Length 


1888 bytes 
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34. SUBROUTINE QINT (DECK SLBC) 
a. Alyorithm 
Perform a quajratic interpolation with the given values. 
b. Input/Output 


None 


' 
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None 

ad, Subroutines Required 
None 

e. Argument List 
{Q1, G2, Q3, Q4) 

f, Length 
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SUBROUTINE CARD (DECK CARD) 


This routine re-ds geometry data from. a tape unit and 
punches the information on cards. 


Alsarithm 

The version for use with the IBM 360 reads the data 
fxom Tape 8 and writes the same information on the 
punch uart (Unit 7). 

Input/Output 


Reads geometry data from Tape 8 and writes the same 
information on the punch unit (Unit 7) 


Error 

None 

Subroutines Required 
None 

Argumen’ List 

None 

Lencth 


469 bytes 
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APPENDIX B 


PROGRAM MNEMONIC LIST 


This appendix contains an alphabetic list of ali symbols used in the Marl: II 
Mod 0 pregram. The list is divided into five fields which are described 
as follows: 


(i) The first field contains the symbol 


(ii) The second field contains the letters I, L, or R, denoting integer, 
logical, or real variable respectively. 





{iii) The third field contains the letters A, C, D, U, denoting argument 
list. common, dimensioned, or urdimensioned, variable respectively. 
The hierarchy of the above letters is A, C, D, U. 


{iv) The fourth field contai.s the definition of the symbol. 


(vy) The fifth field contains the name of the subroutine in which the 
symbol occurs. 
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APPENDIX C 


PROGRAM ARRAYS 


This program contains a number cf subscripted variable arrays. Most of 
these are single parameter arrays and are described in Appendix B. How- . 
ever, some of the arrays are used to represent several different program 
parameters. The rnost important of these are listed below. All of these 
arrays are used in the skin friction calculations. 


Array 
Item Description 
ANGLE(L) Angie through which flow is compressed or expanded. 
For an upper surface, + for expansion 
- for compression 
For alower surface, + for compression 
~ for expansion 
ANGLE(2) Compression or expansion angle (absolute value of ANGLE(1)) 
ANGLE(3) Shock angle for compression and Mach angle for expansion 
TR{1) Altitude 
TR(2) Mach number 
TR(3) Velocity 
TR(4) Angle of attack of flight reference plane 
TR(5) Wall temperature, degrees Rankine, laminar 
TR(6) Wall temperature, degrees Rankine, turbulent 
TR(7) Wall enthalpy, laminar ‘ 
TR(8) Wall enthalpy, turbulent 
TRI{9) Adiabatic wall enthalpy, laminar 
TRI) Adiabatic wall enthalpy, turbulent 
FS(I) Flow conditions before shock ar expansion 
FS(1) Density, slugs/ft? 
FS(2) Pressure, pounds /ft# 
FS(3) Temperature, degrees Rankine 
FS(4) Speed of sound, feet/zec 
FS(5) . Viscosity, slugs/ft-sec | 


C-1 





Array 
Item 


FS(6) 
’ FS(7) 
FS(8) 
BS{I) 


SCF(1) 


SCF(2) 
CFL.(I) 


CFL(1) 


CFL(2) 
CFL(3) 
CFL(4) 


CFL(5) 
CFL(6) 
CFT(D) 
TS(1) 
TS(2) 
RE(l) 
RE(2) 
RF (1) 
RF(2) 
RT(1) 
RT (2) 


Description 
Mach. timber 
Velocity, feet/soc 
Reynolds number per foot 


Flow conditions behind shock or expansion. See FS above 
for individual parameters. 


Total skin friction coefficient based on free stream pro- 
perties and reference area. This is the sum of the proper 
combination of laminar and turbulent coefficients. Value 
in the lift direction. 

Total skin friction coefficient, value in the drag direction. 


Laminar skin friction values. 


Local average coefficient based on incompressible 
relationa, 


Local average coefficient based on local flow conditions. 
Free stream skin friction coefficient based on local length. 


Total skin friction coefficient based on vehicle reference 
area, 


Total:-skin friction coefficient in the lift direction. 

Total skin friction coefficient in the drag direction. 
Turbulent skin friction values (see CFL above). 

Reference temperature, degrees Rankine, laminar. 
Reference temperature, degrees Rankine, turbulent. 
Reference Reynolds number based on local length, laminar. 
Reference Reynolds number based on local length, turbulent. 
Recovery factor, laminar. 

Recovery factor, turbulent. 


Recovery temperature, laminar. 


Recovery temperature, turbulent. 
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project code name, geogrephic location, may be used ao key 
words but will ba followed by an indication of technical con- 
text The aasigoment of links, raises, and weighte ia epuonal. 


Bs 





Unclassified 


mewn 


Security Ciauaification 


a ce ee RR 











